Question about pareto front

[A-ep93]

Hello I am using botorch to optimize my problem that has two outputs and I want to maximize both outputs.
I created my code similar to the tutorial posted here and here are my Pareto front results:

Do you think this Pareto front makes sense? Any advice/ comments?
Thanks in advance.

[saitcakmak]

Hi @A-ep93! That does not look like a Pareto frontier if you're optimizing both objectives. The plot includes many points that are clearly dominated by others. I suspect the plot shows all evaluations, which is what the bottom plots do in the notebook you linked to.

Judging by these plots, and assuming those are the evaluations, it seems that the two objectives are quite strongly correlated here (maybe you don't really need to do multi-objective optimization?). The PF here would be a handful of points that are at upper right corner of your plots. If you want to filter out the Pareto points, you can use the following where Y is the n x m-dim tensor of all evaluations:

from botorch.utils.multi_objective.pareto import is_non_dominated
pareto_mask = is_non_dominated(Y)
pareto_Y = Y[pareto_mask]

[A-ep93]

Thanks for your reply @saitcakmak. Yes, the plot shows all evaluations. Also, I see the points you mentioned about the correlation between outputs and I will surely look into that. But I encountered one other problem:
I changed the dimension of my problem to 4 (it was 1 in previous results I posted). Now my problem has 4 dimensions and two objectives that I want to maximize. Here are the results:

These results are worse than the results I posted previously. Any suggestions to somehow increase the performance?

I post the main parts of my code here (I didn't bring unnecessary parts here) if you want to take a look at what exactly I am doing:

# Here are the bounds of input parameters and the reference point


bounds=torch.tensor([[1,0.1,1000,1e-05],[50,5,30000,5e-04]],dtype=torch.float64)
standard_bounds = torch.zeros(2, 4, **tkwargs)
standard_bounds[1] = 1
ref=torch.tensor([10,20],dtype=torch.float)



def problem(input):
    #This is my black box function
    return torch.tensor(result,dtype=torch.float64)




def generate_initial_data(n=6):
    # generate training data
    train_x = draw_sobol_samples(
        bounds=bounds,n=1, q=n, seed=torch.randint(1000000, (1,)).item()
    ).squeeze(0)
    train_obj = problem(train_x)
    return train_x, train_obj


def initialize_model(train_x, train_obj):
    models = [
        SingleTaskGP(
            train_x,
            train_obj[..., i : i + 1],
            outcome_transform=Standardize(m=1),
        )
        for i in range(train_obj.shape[-1])
    ]
    model = ModelListGP(*models)
    mll = SumMarginalLogLikelihood(model.likelihood, model)
    return mll, model



BATCH_SIZE = 4 if not SMOKE_TEST else 2
NUM_RESTARTS = 20 if not SMOKE_TEST else 2
RAW_SAMPLES = 1024 if not SMOKE_TEST else 4



def optimize_qehvi_and_get_observation(model, train_obj, sampler):
    """Optimizes the qEHVI acquisition function, and returns a new candidate and observation."""
    # partition non-dominated space into disjoint rectangles
    partitioning = NondominatedPartitioning(ref_point=ref, Y=train_obj)
    acq_func = qExpectedHypervolumeImprovement(
        model=model,
        ref_point=ref.tolist(),  # use known reference point 
        partitioning=partitioning,
        sampler=sampler,
    )
    # optimize
    candidates, _ = optimize_acqf(
        acq_function=acq_func,
        bounds=standard_bounds,
        q=BATCH_SIZE,
        num_restarts=NUM_RESTARTS,
        raw_samples=RAW_SAMPLES,  # used for intialization heuristic
        options={"batch_limit": 5, "maxiter": 200, "nonnegative": True},
        sequential=True,
    )
    # observe new values 
    new_x =  unnormalize(candidates.detach(), bounds=bounds)
    new_obj = problem(new_x)
    return new_x, new_obj


def optimize_qparego_and_get_observation(model, train_obj, sampler):
    """Samples a set of random weights for each candidate in the batch, performs sequential greedy optimization 
    of the qParEGO acquisition function, and returns a new candidate and observation."""
    acq_func_list = []
    for _ in range(BATCH_SIZE):
        weights = sample_simplex(2, **tkwargs).squeeze()
        objective = GenericMCObjective(get_chebyshev_scalarization(weights=weights, Y=train_obj))
        acq_func = qExpectedImprovement(  # pyre-ignore: [28]
            model=model,
            objective=objective,
            best_f=objective(train_obj).max(),
            sampler=sampler,
        )
        acq_func_list.append(acq_func)
    # optimize
    candidates, _ = optimize_acqf_list(
        acq_function_list=acq_func_list,
        bounds=standard_bounds,
        num_restarts=NUM_RESTARTS,
        raw_samples=RAW_SAMPLES,  # used for intialization heuristic
        options={"batch_limit": 5, "maxiter": 200},
    )
    # observe new values 
    new_x =  unnormalize(candidates.detach(), bounds=bounds)
    new_obj = problem(new_x)
    return new_x, new_obj



warnings.filterwarnings('ignore', category=BadInitialCandidatesWarning)
warnings.filterwarnings('ignore', category=RuntimeWarning)

N_TRIALS = 3 if not SMOKE_TEST else 2
N_BATCH = 100 if not SMOKE_TEST else 3
MC_SAMPLES = 128  if not SMOKE_TEST else 16

verbose = False

hvs_qparego_all, hvs_qehvi_all, hvs_random_all = [], [], []

hv = Hypervolume(ref_point=ref)


# average over multiple trials
for trial in range(1, N_TRIALS + 1):
    torch.manual_seed(trial)
    
    print(f"\nTrial {trial:>2} of {N_TRIALS} ", end="")
    hvs_qparego, hvs_qehvi, hvs_random = [], [], []
    
    # call helper functions to generate initial training data and initialize model
    train_x_qparego, train_obj_qparego = generate_initial_data(n=6)
    mll_qparego, model_qparego = initialize_model(train_x_qparego, train_obj_qparego)
    
    train_x_qehvi, train_obj_qehvi = train_x_qparego, train_obj_qparego
    train_x_random, train_obj_random = train_x_qparego, train_obj_qparego
    # compute hypervolume 
    mll_qehvi, model_qehvi = initialize_model(train_x_qehvi, train_obj_qehvi)
    
    # compute pareto front
    pareto_mask = is_non_dominated(train_obj_qparego)
    pareto_y = train_obj_qparego[pareto_mask]
    # compute hypervolume
    
    volume = hv.compute(pareto_y)
    
    hvs_qparego.append(volume)
    hvs_qehvi.append(volume)
    hvs_random.append(volume)
    
    # run N_BATCH rounds of BayesOpt after the initial random batch
    for iteration in range(1, N_BATCH + 1):    
        
        t0 = time.time()
        
        # fit the models
        fit_gpytorch_model(mll_qparego)
        fit_gpytorch_model(mll_qehvi)
        
        # define the qEI and qNEI acquisition modules using a QMC sampler
        qparego_sampler = SobolQMCNormalSampler(num_samples=MC_SAMPLES)
        qehvi_sampler = SobolQMCNormalSampler(num_samples=MC_SAMPLES)
        
        # optimize acquisition functions and get new observations
        new_x_qparego, new_obj_qparego = optimize_qparego_and_get_observation(
            model_qparego, train_obj_qparego, qparego_sampler
        )
        new_x_qehvi, new_obj_qehvi = optimize_qehvi_and_get_observation(
            model_qehvi, train_obj_qehvi, qehvi_sampler
        )
        new_x_random, new_obj_random = generate_initial_data(n=BATCH_SIZE)
                
        # update training points
        train_x_qparego = torch.cat([train_x_qparego, new_x_qparego])
        train_obj_qparego = torch.cat([train_obj_qparego, new_obj_qparego])

        train_x_qehvi = torch.cat([train_x_qehvi, new_x_qehvi])
        train_obj_qehvi = torch.cat([train_obj_qehvi, new_obj_qehvi])
    
        train_x_random = torch.cat([train_x_random, new_x_random])
        train_obj_random = torch.cat([train_obj_random, new_obj_random])
        

        # update progress
        for hvs_list, train_obj in zip(
            (hvs_random, hvs_qparego, hvs_qehvi), 
            (train_obj_random, train_obj_qparego, train_obj_qehvi),
        ):
            # compute pareto front
            pareto_mask = is_non_dominated(train_obj)
            pareto_y = train_obj[pareto_mask]
            # compute hypervolume
            volume = hv.compute(pareto_y)
            hvs_list.append(volume)

        # reinitialize the models so they are ready for fitting on next iteration
        # Note: we find improved performance from not warm starting the model hyperparameters
        # using the hyperparameters from the previous iteration
        mll_qparego, model_qparego = initialize_model(train_x_qparego, train_obj_qparego)
        mll_qehvi, model_qehvi = initialize_model(train_x_qehvi, train_obj_qehvi)
        
        t1 = time.time()
        
        if verbose:
            print(
                f"\nBatch {iteration:>2}: Hypervolume (random, qParEGO, qEHVI) = "
                f"({hvs_random[-1]:>4.2f}, {hvs_qparego[-1]:>4.2f}, {hvs_qehvi[-1]:>4.2f}), "
                f"time = {t1-t0:>4.2f}.", end=""
            )
        else:
            print(".", end="")
   
    hvs_qparego_all.append(hvs_qparego)
    hvs_qehvi_all.append(hvs_qehvi)
    hvs_random_all.append(hvs_random)
    
    
    


fig, axes = plt.subplots(1, 3, figsize=(17, 5))
algos = ["Sobol", "qParEGO", "qEHVI"]
cm = plt.cm.get_cmap('viridis')

batch_number = torch.cat(
    [torch.zeros(6), torch.arange(1, N_BATCH+1).repeat(BATCH_SIZE, 1).t().reshape(-1)]
).numpy()
for i, train_obj in enumerate((train_obj_random, train_obj_qparego, train_obj_qehvi)):
    sc = axes[i].scatter(
        train_obj[:, 0].cpu().numpy(), train_obj[:,1].cpu().numpy(), c=batch_number, alpha=0.8,
    )
    axes[i].set_title(algos[i])
    axes[i].set_xlabel("Objective 1")
axes[0].set_ylabel("Objective 2")
norm = plt.Normalize(batch_number.min(), batch_number.max())
sm =  ScalarMappable(norm=norm, cmap=cm)
sm.set_array([])
fig.subplots_adjust(right=0.9)
cbar_ax = fig.add_axes([0.93, 0.15, 0.01, 0.7])
cbar = fig.colorbar(sm, cax=cbar_ax)
cbar.ax.set_title("Iteration")

        

Any suggestions/comments are greatly appreciated.

[saitcakmak]

A few suggestions based on the code:

• Your bounds are quite non-standard. The priors for the SingleTaskGP hyper-parameters may not behave too well in this case. I'd recommend passing in input transform Normalize to your model so that this is nicely handled under the hood.

        SingleTaskGP(
            train_x,
            train_obj[..., i : i + 1],
            outcome_transform=Standardize(m=1),
            input_transform=Normalize(d=4),  # add this line, import from botorch.models.transforms.input
        )

• The reference point you provide is [10, 20] but none of the points in the plots is anywhere near dominating [10, 20]. For NEHVI to work well, you want your samples to have non-zero hyper volume. So, making the reference point achievable would help here. Based on the plots [0, 0] might be a decent starting point.
• There's a number of places where standard bounds are used, then the x is unnormalized to original bounds. These are in general hard to keep track of and there may be a bug hiding behind these. You can change the optimization bounds to your original bounds and remove the unnormalize(candidates, ...) call. If you add the input transform Normalize, this should lead to equivalent behavior, while eliminating a possible place for a bug. Currently, it looks like the optimization is over standard_bounds while the model is defined with inputs in bounds, which may be an explanation for the weird results.

Hope this helps!

[A-ep93]

Thanks for your response @saitcakmak.
I did what you mentioned. Here are the new results:

Do you have any further suggestions?
Also, do you think increasing N_TRIALS, N_BATCH, MC_SAMPLES, RAW_SAMPLES, NUM_RESTARTS, BATCH_SIZE, etc. can increase the performance? If yes, which one do you think is more critical?

I am not sure whether I applied your third comment correctly to my code. So basically, based on your first suggestion, I added input_transform=Normalize(d=4) to my initialize_model function. I then replaced every standard_bounds with bounds in the optimize_acqf functions. Also, I replaced the unnormalize(candidates.detach(), bounds=bounds) comment with candidates.detach(). Did I apply your third comment correctly?

[Balandat]

So on a high level it seems that the outcomes that you're optimizing here are quite correlated. So it's not really clear that multi-objective optimization really provides a lot of value here, since there don't seem to be really any interesting tradeoffs, and any scalarization of the objectives should work here (in which case you could just use a single-objective acquisition function like NEI).

Another meta-comment: Unless you have a need to drop down to the details, I would suggest running the optimization via Ax rather than in BoTorch (see this tutorial). As Sait mentioned above there are a lot of things that can be configured wrong when you are running this in a barebones BoTorch setup - the Ax version takes care of all the transforms and other relevant things on the backend so that you don't have to worry about those details.

[A-ep93]

Thanks for your reply @Balandat. I will surely consider using Ax in my future runs. However, I have one other question now,
as you mentioned, my previous outputs are correlated, so I considered changing them. I chose two other outputs that I believe should not be dependent. Here are my new results:

Now I have two questions:

1. Do you think the correlation problem is solved here?
2. There is no Pareto curve in the outputs. The results seem to be a cloud of points with no specific trends. Any suggestions to improve this behavior to be able to achieve a clear Pareto curve? My code is mostly similar to what I posted previously at the first post of this topic and I had applied the changes suggested by @saitcakmak as well.

[Balandat]

Hi @A-ep93

Do you think the correlation problem is solved here?

I am not sure why it would be a correlation "problem". At the end of the day this is your application / problem and you will know best which outcomes to optimize. If the outcomes you want to optimize are all highly correlated in the right way (e.g. positively if you want to maximize all of them), then that's not a problem - quite the opposite, it's a luxury.

Any suggestions to improve this behavior to be able to achieve a clear Pareto curve?

My response here would be similar to the one above. It almost feels like you are in search of a problem here. The correlations of the outcomes are whatever they are in your problem, and you will have to work with that. Now, regarding the Pareto Frontier: All these plots do have one, it may just not be very "nice" and a natural tradeoff, but there (by definition) is a PF in each of them. You can also see that the qEHVI one has a lot more points with high values for obj1, and these are all points obtained later during the optimzation. So the algorithm seems to be working, it's just that there appears to be a pretty weak dependence of obj 2 on the parameters.

[A-ep93]

Thank you very much for your helpful explanation @Balandat.