Bayesian Quadrature using different kernels

[adrenel1ne]

Hi there!

I am using emukits 0.4.11 version for Bayesian Quadrature! I want to try different kernels and also (if thats possible) lineair combinations of them and also optimize their hyperparameters. I had something like the following:

kernel2 = GPy.kern.src.static.White(2)
kernel3 = GPy.kern.RatQuad(2,  ARD=True, power=1)
kernel4 = GPy.kern.RBF(2, lengthscale=0.2, variance=1)
kernel5 = GPy.kern.Bias(2, variance=100)
kernel = kernel1 + kernel2 + kernel3+ kernel4 + kernel5
gpy_model = GPy.models.GPRegression(X=X_data, Y=force_data, kernel=kernel)
emukit_rbf = RBFGPy(gpy_model.kern)
emukit_measure = LebesgueMeasure.from_bounds(bounds=[(lb_1, ub_1), (lb_2, ub_2)])
emukit_qrbf = QuadratureRBFLebesgueMeasure(emukit_rbf, emukit_measure)
emukit_model = BaseGaussianProcessGPy(kern=emukit_qrbf, gpy_model=gpy_model)
emukit_method = VanillaBayesianQuadrature(base_gp=emukit_model, X=X_data, Y=force_data)

which does run (and also the emukit_method.optimize() does run, however, it does not accept any parameters (I want to set verbose=True).

However, this method uses the RBFGPy to wrap the GPy kernel to an emukit kernel - is this possible? I wasnt so sure, so I also tried:

emukit_measure = LebesgueMeasure.from_bounds(bounds=[(lb_1, ub_1), (lb_2, ub_2)])
emukit_model = create_emukit_model_from_gpy_model(gpy_model, measure=emukit_measure)


emukit_method = VanillaBayesianQuadrature(base_gp=emukit_model, X=X_data, Y=force_data)

,
but here, the create_emukit_model_from_gpy_model cannot handle the summation of the different kernels.

But the first method looks working, but I just wanted to be sure that it is not turning my combination of kernels into another RBF kernel

[adrenel1ne]

I believe I might need to tag @mmahsereci here as well :) , thanks in advance!

[adrenel1ne]

(I am tagging you again, but please enjoy your holiday @mmahsereci )
So here, I am basically doing exactly the same as in my other question, but I want to explore some different kernels as well , so for example, instead of : kernel = GPy.kern.RBF(input_dim=2, lengthscale=0.5, variance=1.0)

I would like to:

kernel1 = GPy.kern.RBF(2, lengthscale=(0.3,0.3), variance=1, ARD=True)
kernel2 = GPy.kern.src.static.White(2,variance = 1)
kernel3 = GPy.kern.RatQuad(2,  ARD=True, power=1)
kernel4 = GPy.kern.RBF(2, lengthscale=(0.4,0.4), variance=1, ARD=True)
kernel5 = GPy.kern.Bias(2, variance=100)
kernel =  kernel1 + kernel2 + kernel3+  kernel4+kernel5

in the second case, can I just use

gpy_model = GPy.models.GPRegression(X=X_data, Y=derivative_data, kernel=kernel)
 emukit_rbf = RBFGPy(gpy_model.kern)
emukit_measure = LebesgueMeasure.from_bounds(bounds=[(lb_1, ub_1), (lb_2, ub_2)])
emukit_qrbf = QuadratureRBFLebesgueMeasure(emukit_rbf, emukit_measure)
emukit_model = BaseGaussianProcessGPy(kern=emukit_qrbf, gpy_model=gpy_model)
emukit_measure = LebesgueMeasure.from_bounds(bounds=[(lb_1, ub_1), (lb_2, ub_2)])
 emukit_model = create_emukit_model_from_gpy_model(gpy_model, measure=emukit_measure)
emukit_method = VanillaBayesianQuadrature(base_gp=emukit_model, X=X_data, Y=derivative_data)

Since I was thinking this is made for an rbf kernel, so I thought maybe I should use emukit_model = create_emukit_model_from_gpy_model(gpy_model, measure=emukit_measure) or something like that, but I was not sure.

Another question that raised during this kernel exploration, was how to use emukit_method.optimize(), and if there is any way to see how it changes the hyperparameters of the kernel.

Thanks in advance!

[adrenel1ne]

Hi again, @mmahsereci ! I have also been working on this, but I do not think I have fixed this. So, I have been using Bayesian Quadrature and I would like to use an other kernel than the rbf kernel. I have seen that you support wrappers for the Matern kernels as well! However, I would also like to use the StandardPeriodic kernel - and also maybe some combinations of kernels (so adding a RBF and a White kernel). I have looked at the documentation, and I think the only thing I need to create is a wrapper. Is this correct? And if so, is something like the following lines correct:

    """
    Wrapper for a sum of GPy RBF and White kernels to be used in Emukit quadrature.

    :param gpy_rbf: An RBF kernel from GPy.
    :param gpy_white: A White kernel from GPy.
    """

    def __init__(self, gpy_kernel):
        
        gpy_rbf = gpy_model.kern.parts[0]
        gpy_white = gpy_model.kern.parts[1]
        self.gpy_rbf = gpy_rbf
        self.gpy_white = gpy_white
        self.gpy_kernel = gpy_rbf + gpy_white  # Sum of kernels

    @property
    def lengthscales(self) -> np.ndarray:
       if self.gpy_rbf.ARD:
           return self.gpy_rbf.lengthscale.values
       return np.full((self.gpy_rbf.input_dim,), self.gpy_rbf.lengthscale[0])

    @property
    def variance(self) -> float:
        """Returns the variance of the total"""
        return self.gpy_rbf.variance.values[0] # I dont think we need to add the variance of the white


    def K(self, x1: np.ndarray, x2: np.ndarray) -> np.ndarray:
        """Computes the full kernel matrix (RBF + White)."""
        return self.gpy_kernel.K(x1, x2)

    def dK_dx1(self, x1: np.ndarray, x2: np.ndarray) -> np.ndarray:
        scaled_vector_diff = np.swapaxes((x1[None, :, :] - x2[:, None, :]) / self.lengthscales**2, 0, -1)
        return -self.K(x1, x2)[None, ...] * scaled_vector_diff


    def dKdiag_dx(self, x: np.ndarray) -> np.ndarray:
        return np.zeros((x.shape[1], x.shape[0]))


# Wrap the GPy model
kernel1 = GPy.kern.RBF(2, lengthscale=2.0, variance=1)
kernel2 = GPy.kern.src.static.White(2,variance = 0.5)

kernel = kernel1 + kernel2

gpy_model = GPy.models.GPRegression(X=X_data, Y=derivative_data, kernel=kernel)
emukit_kernel = SumRBFWhiteGPy(gpy_model.kern)

this seems to work as I expect it to work. However, for the standard periodic, I tried the following:

class IStandardPeriodic(IStandardKernel):
    r"""Interface for a standard periodic kernel.

    .. math::
        k(x, x') = \sigma^2 \exp\left(-\frac{2\sin^2(\pi (x - x') / p)}{\lambda^2} \right),

    where:
    - :math:`\lambda` is the lengthscale,
    - :math:`p` is the period,
    - :math:`\sigma^2` is the variance.

    .. note::
        This is an interface for wrapping standard periodic kernels.
    """

    @property
    def lengthscales(self) -> np.ndarray:
        """Returns the lengthscales of the periodic kernel."""
        raise NotImplementedError

    @property
    def variance(self) -> float:
        """Returns the variance of the periodic kernel."""
        raise NotImplementedError

    @property
    def period(self) -> float:
        """Returns the period of the periodic kernel."""
        raise NotImplementedError

    def dK_dx1(self, x1: np.ndarray, x2: np.ndarray) -> np.ndarray:
        """Gradient of the periodic kernel w.r.t. x1."""
        factor = np.pi * (x1[None, :, :] - x2[:, None, :])  / self.period
        sin_term = np.sin(factor)
        scaled_diff = (4 * np.pi* sin_term * np.cos(factor)) / (self.lengthscales**2 * self.period)
        return -self.K(x1, x2)[None, ...] * np.swapaxes(scaled_diff, 0, -1)

    
    def dKdiag_dx(self, x: np.ndarray) -> np.ndarray:
        """Gradient of the kernel diagonal (variance) w.r.t. x."""
        return np.zeros((x.shape[1], x.shape[0]))  # No change in variance because sin(0) = 0


class wrap_StandardPeriodicGPy(IStandardPeriodic):
    """
    Wrapper for GPy's standard Periodic kernel to be used in Emukit quadrature.

    :param gpy_periodic: A standard periodic kernel from GPy.
    """

    def __init__(self, gpy_periodic: GPy.kern.StdPeriodic):
        assert isinstance(gpy_periodic, GPy.kern.StdPeriodic), "Input must be a GPy StdPeriodic kernel."
        self.gpy_periodic = gpy_periodic

    @property
    def lengthscales(self) -> np.ndarray:
        """Returns the lengthscales of the Periodic kernel."""
        return np.full((self.gpy_periodic.input_dim,), self.gpy_periodic.lengthscale[0])
    @property
    def variance(self) -> float:
        """Returns the variance of the periodic kernel."""
        return self.gpy_periodic.variance.values[0]

    @property
    def period(self) -> float:
        """Returns the period of the periodic kernel."""
        return self.gpy_periodic.period.values[0]

    def K(self, x1: np.ndarray, x2: np.ndarray) -> np.ndarray:
        """Computes the full periodic kernel matrix."""
        return self.gpy_periodic.K(x1, x2)

    def Kdiag(self, x: np.ndarray) -> np.ndarray:
        """Computes the diagonal of the periodic kernel matrix."""
        return self.gpy_periodic.Kdiag(x)




periodic_kernel = GPy.kern.StdPeriodic(2, lengthscale=1.0, variance=1.0, period=10.0)

gpy_model = GPy.models.GPRegression(X=X_data, Y=derivative_data, kernel=periodic_kernel)
emukit_kernel = wrap_StandardPeriodicGPy(gpy_model.kern)

I might have misinterpreted some things, but the acquisition of the points does not really make sense to me using this kernel. This also might be because I am not sure if the derivative def dK_dx1 is correct in this case. So, my main question, is this indeed the way to go to use other kernels for bayesian quadrature? (And maybe you have an idea to fix the periodic kernel as well).

[mmahsereci]

Hi @adrenel1ne , thanks for the updated code snippet, I'll have a look.