Scoring metrics for multidimensional inputs

[DavidSiretMarques]

I'm trying to do a multidimensional input gaussian regressor, and I wanted to check it's performance with the functions in metrics, but when I've tried, it always gives errors (the two tensors should be the same shape).

Reading the code, I've realized that there seems to be no functions for multidimensional inputs. And the tutorial certainly does not mention anything about that.

I'm asking here because I'm new to this library and don't know much, but is there a simple way to do scoring metrics when your inputs are more than 1-dimensional?

If there isn't, then it might be a good proposition for a new feature...

[gpleiss]

Can you post a code example that reproduces the error? The metrics should work for multidimensional inputs.

[DavidSiretMarques]

I'm trying multidimensional input and output, maybe that's the problem...

Let me figure out what I was doing when I asked the question (as it didn't work I started trying other things) and I'll give you some code

[DavidSiretMarques]

Ok, I've figured it out.

I'm using code from this github repo about using MOBO in physics problems. I've adapted the code for the MOBO using only the analysis 1b (as it was the fastest in my machine). It took me some time to get the hang of it, and maybe I'm not using it right (and that might be part of the problem).

But, from what I've been trying I can give you some code to run and see the problem I'm facing:

#Import classical numpy and torch
import torch as to
import numpy as np
#Plotting
import matplotlib.pyplot as plt
#Random # generator
import random
import time

#Optional for making movies
#import atomai as aoi

# Import GP and BoTorch functions
import gpytorch as gpt
from botorch.models import ModelListGP
from botorch.models.gpytorch import GPyTorchModel
from gpytorch.distributions import MultivariateNormal
from gpytorch.kernels import ScaleKernel, RBFKernel
from gpytorch.likelihoods import GaussianLikelihood
from gpytorch.means import ConstantMean
from gpytorch.mlls import ExactMarginalLogLikelihood
from gpytorch.constraints import GreaterThan
from gpytorch.models import ExactGP
from gpytorch.metrics import mean_absolute_error

# Specific for MOBO
from botorch.optim.optimize import optimize_acqf, optimize_acqf_list
from botorch.acquisition.objective import GenericMCObjective
from botorch.utils.multi_objective.scalarization import get_chebyshev_scalarization
from botorch.sampling import SobolQMCNormalSampler
# for qEIHV acquistion function (Posteriori method)
from botorch.acquisition.multi_objective.monte_carlo import qExpectedHypervolumeImprovement
from botorch.utils.multi_objective.box_decompositions.non_dominated import NondominatedPartitioning
# for qParEGO acquisition function (Prioiri method)
from botorch.utils.sampling import sample_simplex #This defines weights to objectives randomly
from botorch.acquisition.monte_carlo import qExpectedImprovement
# Important functions for MOBO- select non-dominated Pareto solutions and compute hypervolume
from botorch.utils.multi_objective.pareto import is_non_dominated
from botorch.utils.multi_objective.hypervolume import Hypervolume

#Import optimizer algorithm to train GP hyperparameters
from torch.optim import SGD
from torch.optim import Adam

#other math functions
from scipy.stats import norm

#########################Gaussian Process########################

class SimpleCustomGP(ExactGP, GPyTorchModel):
    _num_outputs = 1  # to inform GPyTorchModel API

    def __init__(self, train_X, train_Y):
        # squeeze output dim before passing train_Y to ExactGP
        super().__init__(train_X, train_Y.squeeze(-1), GaussianLikelihood())
        ## Here we use constant mean, user can change here with different mean function
        self.mean_module = ConstantMean()
        ## Here we use Radial Basis function kernel, user can change here with different kernel function
        self.covar_module = ScaleKernel(
            base_kernel=RBFKernel(ard_num_dims=train_X.shape[-1]),
        )
        self.to(train_X)  # make sure we're on the right device/dtype

    def forward(self, x):
        mean_x = self.mean_module(x)
        covar_x = self.covar_module(x)
        return MultivariateNormal(mean_x, covar_x)

############## Optimize Hyperparameters of GP##############################
def optimize_hyperparam_trainGP(train_X_norm, train_Y_norm):
    # Gp model fit

    gp_surro1 = SimpleCustomGP(train_X_norm, train_Y_norm[:, 0])
    gp_surro2 = SimpleCustomGP(train_X_norm, train_Y_norm[:, 1])
    gp_surro = ModelListGP(gp_surro1, gp_surro2)
    gp_surro = gp_surro.double()
    gp_surro.likelihood.likelihoods[0].noise_covar.register_constraint("raw_noise", GreaterThan(1e-1))
    gp_surro.likelihood.likelihoods[1].noise_covar.register_constraint("raw_noise", GreaterThan(1e-1))
    mll1 = ExactMarginalLogLikelihood(gp_surro.likelihood.likelihoods[0], gp_surro.models[0])
    mll2 = ExactMarginalLogLikelihood(gp_surro.likelihood.likelihoods[1], gp_surro.models[1])
    # fit_gpytorch_model(mll)
    mll1 = mll1.to(train_X_norm)
    mll2 = mll2.to(train_X_norm)
    gp_surro.models[0].train()
    gp_surro.likelihood.likelihoods[0].train()
    gp_surro.models[1].train()
    gp_surro.likelihood.likelihoods[1].train()
    ## Here we use Adam optimizer with learning rate =0.2, user can change here with different algorithm and/or learning rate for each GP
    optimizer1 = Adam([{'params': gp_surro.models[0].parameters()}], lr=0.2)
    optimizer2 = Adam([{'params': gp_surro.models[1].parameters()}], lr=0.2)

    #How many itereations for finding best hyperparams
    NUM_EPOCHS = 10  
    

    for epoch in range(NUM_EPOCHS):
        # clear gradients
        optimizer1.zero_grad()
        optimizer2.zero_grad()
        # forward pass through the model to obtain the output MultivariateNormal
        output1 = gp_surro.models[0](train_X_norm)
        output2 = gp_surro.models[1](train_X_norm)
        # Compute negative marginal log likelihood
        loss1 = - mll1(output1, gp_surro.models[0].train_targets)
        loss2 = - mll2(output2, gp_surro.models[1].train_targets)
        # back prop gradients
        loss1.backward(retain_graph=True)
        loss2.backward(retain_graph=True)
        # print last iterations
        if (epoch + 1) > NUM_EPOCHS: #Stopping the print for now
            print("GP Model trained:")
            print("Iteration:" + str(epoch + 1))
            print("Loss for ob1:" + str(loss1.item()))
            print("Loss for ob2:" + str(loss2.item()))
            print("noise for obj1:" + str(gp_surro.likelihood.likelihoods[0].noise.item()))
            print("noise for obj2:" + str(gp_surro.likelihood.likelihoods[1].noise.item()))

        optimizer1.step()
        optimizer2.step()

    gp_surro.models[0].eval()
    gp_surro.likelihood.likelihoods[0].eval()
    gp_surro.models[1].eval()
    gp_surro.likelihood.likelihoods[1].eval()
    return gp_surro

# Normalize all data. It is very important to fit GP model and MOBO exploration with normalized data to avoid issues such as
# - decrease of GP performance due to largely spaced real-valued data X, over-weighting on higher valued objectives Y in acquistion function computation.
def normalize_get_initialdata_testfunc(X1, X2, test_X1, test_X2, num_rows, num, problem, testfunction, m):
    # Normalize X
    minX1 = min(X1)
    maxX1 = max(X1)
    minX2 = min(X2)
    maxX2 = max(X2)
    Xmin = np.hstack((minX1, minX2))
    Xmax = np.hstack((maxX1, maxX2))
    X1, X2 = np.meshgrid(X1, X2)
    test_X1, test_X2 = np.meshgrid(test_X1, test_X2)
    X1_norm = (X1 - Xmin[0]) / (Xmax[0] - Xmin[0])
    X2_norm = (X2 - Xmin[1]) / (Xmax[1] - Xmin[1])
    testX1_norm = (test_X1 - Xmin[0]) / (Xmax[0] - Xmin[0])
    testX2_norm = (test_X2 - Xmin[1]) / (Xmax[1] - Xmin[1])

    # Select starting samples randomly as training data
    #reshape to make the points of x1 and x2 a 1D list (instead of 2D list) for easier sampling
    testX1_norm_sampling = np.reshape(testX1_norm, (num_rows ** 2, 1)) 
    testX2_norm_sampling = np.reshape(testX2_norm, (num_rows ** 2, 1))
    test_X1_sampling = np.reshape(test_X1, (num_rows ** 2, 1))
    test_X2_sampling = np.reshape(test_X2, (num_rows ** 2, 1))
    np.random.seed(10)
    ind_start = np.random.choice(len(testX1_norm) - 1, size=num, replace=False)
    train_X = np.zeros((num, 2))
    train_X_norm = np.zeros((num, 2))
    train_y = np.zeros((num, 2))
    t_y = np.zeros((num_rows ** 2, 2))
    m = 0
    print("Function evaluaton for initial randomly selected samples")
    for i in range(0, num):
        train_X_norm[i, 0] = testX1_norm_sampling[ind_start[i]]
        train_X_norm[i, 1] = testX2_norm_sampling[ind_start[i]]
        train_X[i, 0] = test_X1_sampling[ind_start[i]]
        train_X[i, 1] = test_X2_sampling[ind_start[i]]
        print("Function eval #" + str(m + 1))
        if (testfunction == 'func'):
            f1,f2 = arb_func(train_X[i, 0], train_X[i, 1])
        else:
            print("Invalid test function")

        if (problem[0] == 'min'):
            train_y[i, 0] = -f1  # By default, BO settting is for maximization, for minimization we just use f(x) =-f(x)
            t_y[ind_start[i], 0] = -f1  # By default, BO settting is for maximization, for minimization we just use f(x) =-f(x)
        else:
            train_y[i, 0] = f1
            t_y[ind_start[i], 0] = f1

        if (problem[1] == 'min'):
            train_y[i, 1] = -f2  # By default, BO settting is for maximization, for minimization we just use f(x) =-f(x)
            t_y[ind_start[i], 1] = -f2  # By default, BO settting is for maximization, for minimization we just use f(x) =-f(x)
        else:
            train_y[i, 1] = f2
            t_y[ind_start[i], 1] = f2

        # testX1_norm[ind_start[i]] =-1
        # testX2_norm[ind_start[i]] =-1
        m = m + 1
    t_y = np.reshape(t_y, (num_rows, num_rows, 2))

    # Convert to torch
    train_X_norm = to.from_numpy(train_X_norm)
    train_X_norm = train_X_norm.double()
    train_Y = to.from_numpy(train_y)
    
    # Normalize Y (we have 2 Y, if one is much bigger than the other, that could mean trouble)
    testY1_norm = (train_Y[:, 0] - train_Y[:, 0].min()) / (train_Y[:, 0].max() - train_Y[:, 0].min())
    testY2_norm = (train_Y[:, 1] - train_Y[:, 1].min()) / (train_Y[:, 1].max() - train_Y[:, 1].min())
    train_Y_norm = to.empty(train_Y.__len__(), 2)
    train_Y_norm[:, 0] = testY1_norm
    train_Y_norm[:, 1] = testY2_norm

    return X1_norm, X2_norm, testX1_norm, testX2_norm, train_X, train_X_norm, train_Y, train_Y_norm, Xmin, Xmax, t_y, m

##Toy model
def arb_func(x1,x2):
    f1 = x1**2 - 2 * x2 #np.sin(x1)**3-np.cos(x2)**3
    f2 = x1 + x2
    return f1,f2

def MOBO_testfunc(X1, X2, test_X1, test_X2, num_rows, N, BATCH_SIZE, param1, param2, obj1, obj2, acquisition, problem, testfunction):
        ##MOBO framework#######
    m=0
    #normalize all data
    X1_norm, X2_norm, testX1_norm, testX2_norm, train_X, train_X_norm, train_Y, train_Y_norm, Xmin, Xmax, t_y, m=\
        normalize_get_initialdata_testfunc(X1, X2, test_X1, test_X2, num_rows, num, problem, testfunction, m)

    X1, X2 = np.meshgrid(X1, X2)
    test_X1, test_X2 = np.meshgrid(test_X1, test_X2)

    ## Gp model fit
    #Calling function to fit and optimizize Hyperparameter of Gaussian Process (using Adam optimizer)
    # Input args- Torch arrays of normalized training data, parameter X and objective eval Y
    # Output args- Gaussian process model lists
    gp_surro= optimize_hyperparam_trainGP(train_X_norm, train_Y_norm)
    
    mae = mean_absolute_error([gp_surro.models[0](test_X1),gp_surro.models[1](test_X2)],arb_func(test_X1,test_X2))
    print(mae)


#########################################################################################################################################################################################################################################################################
#########################################################################################################################################################################################################################################################################
#########################################################################################################################################################################################################################################################################
# Run MOBO
#########################################################################################################################################################################################################################################################################
#########################################################################################################################################################################################################################################################################
#########################################################################################################################################################################################################################################################################

tic=time.time()

tkwargs = {
    "dtype": to.double,
    "device": to.device("cuda" if to.cuda.is_available() else "cpu"),
}
BATCH_SIZE = 4

# List of options
# 'qEHVI'- Expected Hypervolume acquistion function
# 'qParEGO' -qParEGO
acquisition = 'qEHVI'

#List of control params- Temp, rho, h, deltaG
# Choose any two control params to create the parameter space considering the same order-
param1 = 'X1'
param2 = 'X2'
# Choose any two objectives to do MOO (MOBO)-
obj1 = 'f1'
obj2 = 'f2'

# Test function
testfunction = 'func'

# Input choice of problem- where it is minimization or maximization for objective 1 and 2 respectively
# List of options
# 'min'- Minimization
# 'max' -Maximization (default setting)
problem =['max','max']

#Define parameter space based on above choice of params (2D only)
num_rows =50
num = 10  # Starting samples
#N = (0.1 * num_rows ** 2) - num  # Maximum evaluations
#N = int(N)
N=10

if (testfunction == 'func'):
    X1 = np.linspace(-3, 3, num_rows)
    test_X1 = np.linspace(-3, 3, num_rows)
    X2 = X1
    test_X2 = test_X1
else:
    print("Invalid test function")

#Call functions to run MOBO (classical)
#input args - X1, X2- grid numpy arrays for param X1, X2; test_X1, test_X2- grid numpy arrays for param X1, X2 (for optional posterior calculations); num_rows, N, BATCH_SIZE- explained in text sections; 
#param1, param2 - string names of parameter 1, 2; obj1, obj2 - string names of objective 1, 2; acquisition - string name as provided in text section for user choice of acquistion function method;
# problem - string array as provided in text section for user choice of maxmimization or minimization of objectives; testfunction - string name as provided in text section for user selection of test problems
MOBO_testfunc(X1, X2, test_X1, test_X2, num_rows, N, BATCH_SIZE, param1, param2, obj1, obj2, acquisition, problem, testfunction)

print(f'Time elapsed: {(time.time()-tic)} s.')

It's possible (even probable) that the problem has nothing to do with the scoring function and a lot to do with something I'm not understanding correctly from the code and its usage...

[DavidSiretMarques]

Update: The problem is not with multidimensional inputs, but with multidimensional outputs!!!!

I've come back to this problem after some time, and I've run the code for a single output bayesian optimization without any kind of error, but when I try the Multiple objective Bayesian optimization, I start getting errors...

It might be related with the fact that for MOBO the output is not a MultivariateNormal, but a list of multivariate normals, and that might mess up everything...

[DavidSiretMarques]

Ok, I managed to solve my problem (scoring multidimensional outputs. handling multidimensional inputs is done correctly).

Doing:

    mean_squared_error(model(test_x,test_x)[0], test_y[:,0])
    mean_squared_error(model(test_x,test_x)[1], test_y[:,1])

seems to deal correctly with scoring each output