Multi-objective for Rank-Weighted GP Ensemble (RGPE) Model - Problems with the forward() instantiation of the custom model (MultitaskMultivariateNormal)

[MoBurmeister]

Hello together,

I am trying to setup a RGPE Model for multi-objective problems but I am running into problems when setting up the forward method of the costum RGPE model. I took the RGPE botorch tutorial as an example, but I am still getting stuck when trying to setup the MultitaskMultivariateNormal.

The weights are calculated in a different method are not part of my problem here (right now I have done it via pareto dominance comparison, but I might switch it to obective-wise comparison"). This is my setup:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import torch
import time  
import botorch
from botorch.models import MultiTaskGP, SingleTaskGP, ModelList
from botorch.models.transforms.outcome import OutcomeTransform, Standardize
from botorch.models.transforms.input import InputTransform, Normalize
from gpytorch.mlls.exact_marginal_log_likelihood import ExactMarginalLogLikelihood
from botorch.fit import fit_gpytorch_model, fit_gpytorch_mll
from botorch.acquisition.multi_objective import qNoisyExpectedHypervolumeImprovement
from botorch.optim import optimize_acqf
from gpytorch.mlls import SumMarginalLogLikelihood
from gpytorch.mlls.exact_marginal_log_likelihood import ExactMarginalLogLikelihood
from botorch.models import SingleTaskGP
from botorch.fit import fit_gpytorch_mll
from botorch.models.transforms.outcome import Standardize
from botorch.models.transforms.input import Normalize
import torch
from botorch.sampling.normal import SobolQMCNormalSampler
from botorch.models.gpytorch import GPyTorchModel
from gpytorch.models import GP
from gpytorch.distributions import MultivariateNormal, MultitaskMultivariateNormal
from gpytorch.lazy import PsdSumLazyTensor
from gpytorch.likelihoods import LikelihoodList
from torch.nn import ModuleList
from botorch.acquisition.multi_objective.logei import qLogNoisyExpectedHypervolumeImprovement
from linear_operator.operators import BlockDiagLinearOperator, DenseLinearOperator

#data generation: 


#data generation: 
 
from botorch.test_functions.multi_objective import VehicleSafety

bounds = torch.tensor([[1.0, 1.0, 1.0, 1.0, 1.0], [3.0, 3.0, 3.0, 3.0, 3.0]])

train_X = torch.rand(10, 5, dtype=torch.float64) * (bounds[1] - bounds[0]) + bounds[0]
function = VehicleSafety()


def evaluate_train_y(x):
    # Evaluate the train_Y values for the given x tensor
    y = function.evaluate_true(x)
    
    y = -y
    
    return y

train_Y = evaluate_train_y(train_X)   

print(train_X, train_Y)

train_X = train_X.to(dtype=torch.float64)
train_Y = train_Y*-1
train_Y = train_Y.to(dtype=torch.float64)

model = SingleTaskGP(train_X=train_X, train_Y=train_Y, input_transform=Normalize(d=5), outcome_transform=Standardize(m=3))

mll = ExactMarginalLogLikelihood(likelihood=model.likelihood, model=model) 

mll = fit_gpytorch_model(mll)

model_2 = SingleTaskGP(train_X=train_X, train_Y=train_Y*0.7, input_transform=Normalize(d=5), outcome_transform=Standardize(m=3))

mll_2 = ExactMarginalLogLikelihood(likelihood=model_2.likelihood, model=model_2) 

mll_2 = fit_gpytorch_model(mll_2)

posterior = model.posterior(train_X)

class RGPE(GP, GPyTorchModel):

    def __init__(self, models, weights, num_outputs):
        super().__init__()
        self.models = ModuleList(models)
        for m in models:
            if not hasattr(m, "likelihood"):
                raise ValueError(
                    "RGPE currently only supports models that have a likelihood (e.g. ExactGPs)"
                )
        self.likelihood = LikelihoodList(*[m.likelihood for m in models])
        self.weights = weights
        self._num_outputs = num_outputs
        self.to(weights)

    def forward(self, x):
        # Compute the weighted mean and covariance
        weighted_means = []
        weighted_covariances = []
        
        for i, model in enumerate(self.models):
            posterior = model.posterior(x)
            mean_x = posterior.mean
            covar_x = posterior.mvn.lazy_covariance_matrix
            
            # Compute weighted mean and covariance
            weighted_means.append(self.weights[i] * mean_x)
            weighted_covariances.append(self.weights[i]**2 * covar_x)
        
        # Sum the weighted means and covariances
        mean_x = sum(weighted_means)
        covar_x = PsdSumLazyTensor(*weighted_covariances)
        # Return the MultivariateNormal with the weighted mean and covariance
        return MultitaskMultivariateNormal(mean_x, covar_x)


rgpe = RGPE(models=[model, model_2], weights=torch.tensor([0.5, 0.5], dtype=torch.float64), num_outputs=3)

rgpe_posterior = rgpe.posterior(train_X)



ref_point = torch.tensor([-1864.72022, -11.81993945, -0.2903999384], dtype=torch.float64)  # Two objectives


print("X_baseline shape:", train_X.shape)
print("train_Y shape:", train_Y.shape)
print("train_X shape:", train_X.shape)
print("ref_point shape:", ref_point.shape)
print("model output shape:", rgpe.posterior(train_X).mean.shape) 
print(f"model_1 output shape: {model.posterior(train_X).mean.shape}")

acq_function = qLogNoisyExpectedHypervolumeImprovement(model=rgpe, 
                                                            ref_point=ref_point, 
                                                            X_baseline=train_X,
                                                            prune_baseline=False)
        

candidate, acq_value = optimize_acqf(
acq_function=acq_function,
bounds=bounds,
q=1,
num_restarts=40,
raw_samples=512
        )

print(candidate, acq_value)

I am stuck getting this error:

Traceback (most recent call last):
File "c:\Users\Moritz\Documents\BasicGP_RS\smart_doe_bayesian_optimization\trial_area\testtest.py", line 128, in
acq_function = qLogNoisyExpectedHypervolumeImprovement(model=rgpe,
File "C:\Users\Moritz\anaconda3\envs\TL_GP\lib\site-packages\botorch\acquisition\multi_objective\logei.py", line 420, in init
NoisyExpectedHypervolumeMixin.init(
File "C:\Users\Moritz\anaconda3\envs\TL_GP\lib\site-packages\botorch\utils\multi_objective\hypervolume.py", line 644, in init
self._set_cell_bounds(num_new_points=X_baseline.shape[0])
File "C:\Users\Moritz\anaconda3\envs\TL_GP\lib\site-packages\botorch\utils\multi_objective\hypervolume.py", line 708, in _set_cell_bounds
self._baseline_L = self._compute_root_decomposition(posterior=posterior)
File "C:\Users\Moritz\anaconda3\envs\TL_GP\lib\site-packages\botorch\acquisition\cached_cholesky.py", line 118, in _compute_root_decomposition
lazy_covar = extract_batch_covar(posterior.distribution)
File "C:\Users\Moritz\anaconda3\envs\TL_GP\lib\site-packages\botorch\utils\low_rank.py", line 37, in extract_batch_covar
raise BotorchError(
botorch.exceptions.errors.BotorchError: Expected BlockDiagLinearOperator, but got <class 'linear_operator.operators.psd_sum_linear_operator.PsdSumLinearOperator'>.

Somehow I cannot figure out how to "sum" the different lazy_covariance_matrixes.
When I just "pass" a lazy_covariance_matrix of one of the singletaskgp's, everything seems to work.
I hope you can pinpoint me into the right direction on how to solve the error.

Also I would appreciate some intution about the "unstandardization" that is done in the botorch tutorial:

# unstandardize predictions
posterior_mean = posterior.mean.squeeze(-1) * model.Y_std + model.Y_mean
posterior_cov = posterior.mvn.lazy_covariance_matrix * model.Y_std.pow(2)

What would be an efficient way of integrating it into my application and do I even need it?

Thanks in advance for any help and tips!

[MoBurmeister]

Follow up to my post yesterday:

I figured out a way that seems to work for the addition of the weighted lazy_covariance_matrixes:

def forward(self, x):
        # Compute the weighted mean and covariance
        weighted_means = []
        weighted_covariances = []
        
        for i, model in enumerate(self.models):
            posterior = model.posterior(x)
            mean_x = posterior.mean
            covar_x = posterior.mvn.lazy_covariance_matrix

            #print(f"covar_x {i} dense: {covar_x.to_dense()}\n")
            
            # Compute weighted mean and covariance
            weighted_means.append(self.weights[i] * mean_x)
            weighted_covariances.append(self.weights[i]**2 * covar_x)

            #print(f"Weighted Covariance {i} (Dense): \n{weighted_covariances[i].to_dense()}\n")

        #Sum the weighted means and covariances
        mean_x = sum(weighted_means)
        covar_x = weighted_covariances[0]  # Start with the first covariance
    
        #Sum the weighted covariances - all are BlockDiagLinearOperator
        base_covariances = [covar.base_linear_op for covar in weighted_covariances]

        summed_covariance = BlockDiagLinearOperator(sum(base_covariances))

        #print(f"Summed Covariance (Dense): \n{summed_covariance.to_dense()}\n")


        #Return MultivariateNormal with weighted mean and covariance
        return MultitaskMultivariateNormal(mean_x, summed_covariance)

Now I am able to apply the acquisition function and the optimize_acqf without running into errors.

From the printout this seems to work:

covar_x 0 dense: tensor([[2.0174e-02, 7.5638e-07, 4.3004e-06, 0.0000e+00, 0.0000e+00, 0.0000e+00,
0.0000e+00, 0.0000e+00, 0.0000e+00],
[7.5638e-07, 2.0174e-02, 1.1110e-06, 0.0000e+00, 0.0000e+00, 0.0000e+00,
0.0000e+00, 0.0000e+00, 0.0000e+00],
[4.3004e-06, 1.1110e-06, 2.0174e-02, 0.0000e+00, 0.0000e+00, 0.0000e+00,
0.0000e+00, 0.0000e+00, 0.0000e+00],
[0.0000e+00, 0.0000e+00, 0.0000e+00, 8.2128e-02, 3.0843e-06, 1.7545e-05,
0.0000e+00, 0.0000e+00, 0.0000e+00],
[0.0000e+00, 0.0000e+00, 0.0000e+00, 3.0843e-06, 8.2128e-02, 4.5249e-06,
0.0000e+00, 0.0000e+00, 0.0000e+00],
[0.0000e+00, 0.0000e+00, 0.0000e+00, 1.7545e-05, 4.5249e-06, 8.2128e-02,
0.0000e+00, 0.0000e+00, 0.0000e+00],
[0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00,
2.1836e-05, 8.1880e-10, 4.6666e-09],
[0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00,
8.1880e-10, 2.1836e-05, 1.2055e-09],
[0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00,
4.6666e-09, 1.2055e-09, 2.1836e-05]], dtype=torch.float64)

Weighted Covariance 0 (Dense):
tensor([[5.0434e-03, 1.8909e-07, 1.0751e-06, 0.0000e+00, 0.0000e+00, 0.0000e+00,
0.0000e+00, 0.0000e+00, 0.0000e+00],
[1.8909e-07, 5.0434e-03, 2.7774e-07, 0.0000e+00, 0.0000e+00, 0.0000e+00,
0.0000e+00, 0.0000e+00, 0.0000e+00],
[1.0751e-06, 2.7774e-07, 5.0434e-03, 0.0000e+00, 0.0000e+00, 0.0000e+00,
0.0000e+00, 0.0000e+00, 0.0000e+00],
[0.0000e+00, 0.0000e+00, 0.0000e+00, 2.0532e-02, 7.7107e-07, 4.3862e-06,
0.0000e+00, 0.0000e+00, 0.0000e+00],
[0.0000e+00, 0.0000e+00, 0.0000e+00, 7.7107e-07, 2.0532e-02, 1.1312e-06,
0.0000e+00, 0.0000e+00, 0.0000e+00],
[0.0000e+00, 0.0000e+00, 0.0000e+00, 4.3862e-06, 1.1312e-06, 2.0532e-02,
0.0000e+00, 0.0000e+00, 0.0000e+00],
[0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00,
5.4591e-06, 2.0470e-10, 1.1666e-09],
[0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00,
2.0470e-10, 5.4591e-06, 3.0137e-10],
[0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00,
1.1666e-09, 3.0137e-10, 5.4591e-06]], dtype=torch.float64)

covar_x 1 dense: tensor([[9.8850e-03, 3.7062e-07, 2.1072e-06, 0.0000e+00, 0.0000e+00, 0.0000e+00,
0.0000e+00, 0.0000e+00, 0.0000e+00],
[3.7062e-07, 9.8850e-03, 5.4438e-07, 0.0000e+00, 0.0000e+00, 0.0000e+00,
0.0000e+00, 0.0000e+00, 0.0000e+00],
[2.1072e-06, 5.4438e-07, 9.8850e-03, 0.0000e+00, 0.0000e+00, 0.0000e+00,
0.0000e+00, 0.0000e+00, 0.0000e+00],
[0.0000e+00, 0.0000e+00, 0.0000e+00, 4.0243e-02, 1.5113e-06, 8.5969e-06,
0.0000e+00, 0.0000e+00, 0.0000e+00],
[0.0000e+00, 0.0000e+00, 0.0000e+00, 1.5113e-06, 4.0243e-02, 2.2172e-06,
0.0000e+00, 0.0000e+00, 0.0000e+00],
[0.0000e+00, 0.0000e+00, 0.0000e+00, 8.5969e-06, 2.2172e-06, 4.0243e-02,
0.0000e+00, 0.0000e+00, 0.0000e+00],
[0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00,
1.0700e-05, 4.0121e-10, 2.2866e-09],
[0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00,
4.0121e-10, 1.0700e-05, 5.9069e-10],
[0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00,
2.2866e-09, 5.9069e-10, 1.0700e-05]], dtype=torch.float64)

Weighted Covariance 1 (Dense):
tensor([[2.4713e-03, 9.2656e-08, 5.2679e-07, 0.0000e+00, 0.0000e+00, 0.0000e+00,
0.0000e+00, 0.0000e+00, 0.0000e+00],
[9.2656e-08, 2.4713e-03, 1.3609e-07, 0.0000e+00, 0.0000e+00, 0.0000e+00,
0.0000e+00, 0.0000e+00, 0.0000e+00],
[5.2679e-07, 1.3609e-07, 2.4713e-03, 0.0000e+00, 0.0000e+00, 0.0000e+00,
0.0000e+00, 0.0000e+00, 0.0000e+00],
[0.0000e+00, 0.0000e+00, 0.0000e+00, 1.0061e-02, 3.7782e-07, 2.1492e-06,
0.0000e+00, 0.0000e+00, 0.0000e+00],
[0.0000e+00, 0.0000e+00, 0.0000e+00, 3.7782e-07, 1.0061e-02, 5.5430e-07,
0.0000e+00, 0.0000e+00, 0.0000e+00],
[0.0000e+00, 0.0000e+00, 0.0000e+00, 2.1492e-06, 5.5430e-07, 1.0061e-02,
0.0000e+00, 0.0000e+00, 0.0000e+00],
[0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00,
2.6750e-06, 1.0030e-10, 5.7165e-10],
[0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00,
1.0030e-10, 2.6750e-06, 1.4767e-10],
[0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00,
5.7165e-10, 1.4767e-10, 2.6750e-06]], dtype=torch.float64)

Summed Covariance (Dense):
tensor([[7.5146e-03, 2.8175e-07, 1.6019e-06, 0.0000e+00, 0.0000e+00, 0.0000e+00,
0.0000e+00, 0.0000e+00, 0.0000e+00],
[2.8175e-07, 7.5146e-03, 4.1384e-07, 0.0000e+00, 0.0000e+00, 0.0000e+00,
0.0000e+00, 0.0000e+00, 0.0000e+00],
[1.6019e-06, 4.1384e-07, 7.5146e-03, 0.0000e+00, 0.0000e+00, 0.0000e+00,
0.0000e+00, 0.0000e+00, 0.0000e+00],
[0.0000e+00, 0.0000e+00, 0.0000e+00, 3.0593e-02, 1.1489e-06, 6.5354e-06,
0.0000e+00, 0.0000e+00, 0.0000e+00],
[0.0000e+00, 0.0000e+00, 0.0000e+00, 1.1489e-06, 3.0593e-02, 1.6855e-06,
0.0000e+00, 0.0000e+00, 0.0000e+00],
[0.0000e+00, 0.0000e+00, 0.0000e+00, 6.5354e-06, 1.6855e-06, 3.0593e-02,
0.0000e+00, 0.0000e+00, 0.0000e+00],
[0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00,
8.1341e-06, 3.0500e-10, 1.7383e-09],
[0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00,
3.0500e-10, 8.1341e-06, 4.4905e-10],
[0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00,
1.7383e-09, 4.4905e-10, 8.1341e-06]], dtype=torch.float64)

Any intuition, if I am going into the right direction would be appreciated :)

Also I am still uncertain about the unstandardization. This is still an unclear topic for me.

[sdaulton]

Hi @MoBurmeister!

Since you are treating each outcome (the values in train_Y) independently, a simpler approach might be to construct a RGPE for each outcome (across all tasks --there are 2 tasks in this case). You could use the same weights if desired (or not) and then pass all RGPE models into a ModelListGP to obtain a multi-output GP that outputs the 3 outcomes. This avoids having to handle treating the output/outcome dimension as a batch dimension in the RGPE model. The

Regarding standardization and normalization, those look good to me here. You may want to pass bounds to Normalize since the bounds are known here. The unstandardization is applied to each submodel in the RGPE when model.posterior is called in RGPE.forward.

Aside the RGPE tutorial looks pretty date (calling FixedNoiseGP, which is deprecated), so I'll put up a PR to update that.

[MoBurmeister]

Thanks a lot for your suggestion, the extra step in combining the RGPE into a single ModelList is a good idea and I got it running:

import torch
from botorch.models import SingleTaskGP, ModelListGP
from botorch.models.transforms.outcome import Standardize
from botorch.models.transforms.input import Normalize
from gpytorch.mlls.exact_marginal_log_likelihood import ExactMarginalLogLikelihood
from botorch.fit import fit_gpytorch_model
from botorch.optim import optimize_acqf
from gpytorch.mlls.exact_marginal_log_likelihood import ExactMarginalLogLikelihood
from botorch.models import SingleTaskGP
from botorch.models.transforms.outcome import Standardize
from botorch.models.transforms.input import Normalize
import torch
from botorch.models.gpytorch import GPyTorchModel
from gpytorch.models import GP
from gpytorch.distributions import MultivariateNormal
from gpytorch.lazy import PsdSumLazyTensor
from gpytorch.likelihoods import LikelihoodList
from torch.nn import ModuleList
from botorch.acquisition.multi_objective.logei import qLogNoisyExpectedHypervolumeImprovement

def generate_y(train_X, variation_factor):
    """
    Generate output y based on the input train_X with a variation factor.

    Args:
        train_X (Tensor): Input tensor of shape (n, d) where n is the number of data points
                          and d is the dimension of each data point (here, d=2).
        variation_factor (float): A factor that introduces variation into the output.

    Returns:
        Tensor: Output tensor y of shape (n, 1).
    """
    x1 = train_X[:, 0]
    x2 = train_X[:, 1]
    
    # Sinusoidal function with variation factor applied
    y = torch.sin(x1 * (2 * torch.pi) * variation_factor) + torch.cos(x2 * (2 * torch.pi)) * (1 + variation_factor)
    
    return y.unsqueeze(-1)

# Generate sample two-dimensional training data
train_X = torch.rand(10, 2, dtype=torch.float64)  # 10 training points, each with 2 dimensions

# Variation factors for different models
variation_factors = [1.0, 2.0, 0.5]

# Generate train_Y using the single function with different variation factors
train_Y_1 = generate_y(train_X, variation_factors[0])
train_Y_2 = generate_y(train_X, variation_factors[1])
train_Y_3 = generate_y(train_X, variation_factors[2])

input_transform = Normalize(d=2)

output_transform_1 = Standardize(m=1)

# Initialize the SingleTaskGP models
model_1 = SingleTaskGP(train_X, train_Y_1, input_transform=input_transform, outcome_transform=output_transform_1)
model_2 = SingleTaskGP(train_X, train_Y_2, input_transform=input_transform, outcome_transform=output_transform_1)
model_3 = SingleTaskGP(train_X, train_Y_3, input_transform=input_transform, outcome_transform=output_transform_1)

# Fit the models
mll_1 = ExactMarginalLogLikelihood(model_1.likelihood, model_1)
mll_2 = ExactMarginalLogLikelihood(model_2.likelihood, model_2)
mll_3 = ExactMarginalLogLikelihood(model_3.likelihood, model_3)

fit_gpytorch_model(mll_1)
fit_gpytorch_model(mll_2)
fit_gpytorch_model(mll_3)

# Define bounds for the input space
bounds = torch.tensor([[0.0, 0.0], [1.0, 1.0]], dtype=torch.float64)  

modellist = [model_1, model_2, model_3]

weights = torch.tensor([0.5, 0.3, 0.2], dtype=torch.float64)  # Weights for the models

class RGPE(GP, GPyTorchModel):
    """
    Rank-weighted GP ensemble. Note: this class inherits from GPyTorchModel which provides an
        interface for GPyTorch models in botorch.
    """

    _num_outputs = 1  # metadata for botorch

    def __init__(self, models, weights):
        super().__init__()
        self.models = ModuleList(models)
        for m in models:
            if not hasattr(m, "likelihood"):
                raise ValueError(
                    "RGPE currently only supports models that have a likelihood (e.g. ExactGPs)"
                )
        self.likelihood = LikelihoodList(*[m.likelihood for m in models])
        self.weights = weights
        self.to(weights)

    def forward(self, x):
        weighted_means = []
        weighted_covars = []
        # filter model with zero weights
        # weights on covariance matrices are weight**2
        non_zero_weight_indices = (self.weights**2 > 0).nonzero()
        non_zero_weights = self.weights[non_zero_weight_indices]
        # re-normalize
        non_zero_weights /= non_zero_weights.sum()

        for non_zero_weight_idx in range(non_zero_weight_indices.shape[0]):
            raw_idx = non_zero_weight_indices[non_zero_weight_idx].item()
            model = self.models[raw_idx]
            posterior = model.posterior(x)
            # unstandardize predictions
            posterior_mean = posterior.mean.squeeze(-1) * model.outcome_transform.stdvs + model.outcome_transform.means
            posterior_cov = posterior.mvn.lazy_covariance_matrix * model.outcome_transform.stdvs.pow(2)
            # apply weight
            weight = non_zero_weights[non_zero_weight_idx]
            weighted_means.append(weight * posterior_mean)
            weighted_covars.append(posterior_cov * weight**2)
        # set mean and covariance to be the rank-weighted sum the means and covariances of the
        # base models and target model
        mean_x = torch.stack(weighted_means).sum(dim=0)
        covar_x = PsdSumLazyTensor(*weighted_covars)
        return MultivariateNormal(mean_x, covar_x)

# This part here would be adjusted to catch just the SingleTaskGP models for the considered objective each
rgpe = RGPE(modellist, weights)

rgpe_2 = RGPE(modellist, weights)

list_of_models = [rgpe, rgpe_2]

modellist = ModelListGP(*list_of_models)

# Define the Noisy Expected Hypervolume Improvement (NEHVI) acquisition function

acq_function = qLogNoisyExpectedHypervolumeImprovement(model=modellist, 
                                                       ref_point=torch.tensor([0.0, 0.0], dtype=torch.float64),
                                                       X_baseline=train_X,
                                                       prune_baseline=False)

candidate, acq_value = optimize_acqf(
            acq_function=acq_function,
            bounds=bounds,
            q=1,
            num_restarts=40,
            raw_samples=512
        )


print(f"Optimal solution: {candidate}")

As mentioned in the comment this code is still missing the part of collecting the SingleTaskGPs for each objective and feeding them into the ccorresponding RGPE model. But nevertheless the code here works for the scenario and I will adjust the rest in my on repository.

I changed the unstandardization so that it works with SingleTaskGP:

 # unstandardize predictions
            posterior_mean = posterior.mean.squeeze(-1) * model.outcome_transform.stdvs + model.outcome_transform.means
            posterior_cov = posterior.mvn.lazy_covariance_matrix * model.outcome_transform.stdvs.pow(2)

Hopefully this is doing the right job.

Do you maybe have an intuition on the weight calculation? I am a bit uncertain about the best choice here.
Either I could calculate the weights via pareto dominance and would keep them consistant for all RGPE's. Or I though about calculating the weights individually for each RGPE like it is done in the original approach. Maybe someone has an intuition about what might be the best choice here.

Also there is no way in activating the "prune_baseline" for this setup right?

Thanks a lot for all your help!

[sdaulton]

I changed the unstandardization so that it works with SingleTaskGP:

model.posterior automatically applies the OutputTransform, so this shouldn't be necessary since Standardize is used on the SingleTaskGPs. I merged in a PR yesterday to update the RGPE tutorial to do this.

Do you maybe have an intuition on the weight calculation? I am a bit uncertain about the best choice here.
Either I could calculate the weights via pareto dominance and would keep them consistant for all RGPE's. Or I though about calculating the weights individually for each RGPE like it is done in the original approach. Maybe someone has an intuition about what might be the best choice here.

I would use the original approach. In the RGPE, we want to weight each model according to how well it is able to predict the target task. Doing this for each objective independently is quite intuitive because the optimal weighting (across tasks) may vary across objectives

[sdaulton]

Also there is no way in activating the "prune_baseline" for this setup right?

Prune baseline can still be activated here (by passing prune_baseline=True. This will prune what points are considered as the pareto optimal incumbents in qLogNEHVI.