Own Covariance function for the Heat Equation with pyGPsΒΆ

The definition of this kernel was an attempt to train a pyGPs-GPR-model with the kernel given by \(k_{ff}\) in Chapter 3.1.2, in accordance with the approach as it was described in Chapter 5.

The coefficients were calculated using the heat_equation_with_pygps-v2.ipynb notebook on our GitHub-Page.

The optimal parameters for \(k_{uu}\) were estimated as \([\sigma_u, l_u] = [2e-05, 1291919.81]\).

In [1]:
# The kernel has to be added to cov.py, which is located (with Windows and
# Anaconda) in C:\Users\SurfaceAdmin\Anaconda3\Lib\site-packages\pyGPs\Core

# Unfortunately the kernel is basically equal to zero and the optimizer
# will stick with the default value for the hyperparameter, which in turn
# returns 1.0 for phi.

# Own Kernel function (the try-except block was only added to prevent
# error-messages concerning the unknown Kernel-class after execution of
# the code):

try:
    class MyKernel2(Kernel):

        def __init__(self, log_phi=0.):
            self.hyp = [log_phi]

        def getCovMatrix(self,x=None,z=None,mode=None):
            self.checkInputGetCovMatrix(x,z,mode)
            p = np.exp(self.hyp[0])           # phi
            A = 0
            if not x is None:
                n, D = x.shape
                A = np.zeros((n,n))
            if not z is None:
                nn, D = z.shape
                A = np.zeros((nn,nn))
            if mode == 'self_test':
                A = np.zeros((nn,1))
            elif mode == 'train':             # compute covariance matrix for
                                              # the training set
                A = np.zeros((n,n))
                for i in range(n):
                    for j in range(n):
                        A[i][j] = (p**2*(1.43e-34*(x[i][1] - x[j][1])**4 - \
                                1.11e-27*(x[i][1] - x[j][1])**2 + 7.17e-22) - \
                                2.39e-22*(x[i][0] - x[j][0])**2 + 3.09e-16)* \
                                np.exp(-3.87e-7*(x[i][0] - x[j][0])**2 - \
                                3.87e-7*(x[i][1] - x[j][1])**2)

            elif mode == 'cross':   # compute covariance between data sets x and z
                m = z.shape[0]
                A = np.zeros((n,m))
                for i in range(n):
                    for j in range(m):
                        A[i][j] = (p**2*(1.43e-34*(x[i][1] - z[j][1])**4 - \
                                1.11e-27*(x[i][1] - z[j][1])**2 + 7.17e-22) - \
                                2.39e-22*(x[i][0] - z[j][0])**2 + 3.09e-16)* \
                                np.exp(-3.87e-7*(x[i][0] - z[j][0])**2 - \
                                3.87e-7*(x[i][1] - z[j][1])**2)
            return A

    # We are taking the derivative w.r.t. p, but are multiplying it with 2*p or p,
    # since that seems to be the pattern in the source code of pyGPs as well:

        def getDerMatrix(self,x=None,z=None,mode=None,der=None):
            self.checkInputGetCovMatrix(x,z,mode)
            p = np.exp(self.hyp[0])           # phi
            n = 0
            if not x is None:
                n, D = x.shape
            if not z is None:
                nn, D = z.shape
            if mode == 'self_test':           # self covariances for the test cases
                A = np.zeros((nn,1))
            elif mode == 'train':             # compute covariance matrix for
                                              # the dataset x
                A = np.zeros((n,n))
                for i in range(n):
                    for j in range(n):
                        A[i][j] = (p**2*(1.43e-34*(x[i][1] - x[j][1])**4 - \
                                1.11e-27*(x[i][1] - x[j][1])**2 + 7.17e-22) - \
                                2.39e-22*(x[i][0] - x[j][0])**2 + 3.09e-16)* \
                                np.exp(-3.87e-7*(x[i][0] - x[j][0])**2 - \
                                3.87e-7*(x[i][1] - x[j][1])**2)

            elif mode == 'cross':             # compute covariance between data
                                              # sets x and z
                A = np.zeros((n,nn))
                for i in range(n):
                    for j in range(nn):
                        A[i][j] = (p**2*(1.43e-34*(x[i][1] - z[j][1])**4 - \
                                1.11e-27*(x[i][1] - z[j][1])**2 + 7.17e-22) - \
                                2.39e-22*(x[i][0] - z[j][0])**2 + 3.09e-16)* \
                                np.exp(-3.87e-7*(x[i][0] - z[j][0])**2 - \
                                3.87e-7*(x[i][1] - z[j][1])**2)
            if der == 0:    # compute derivative matrix wrt 1st parameter
                if mode == 'train':
                    A = np.zeros((n,n))
                    for i in range(n):
                        for j in range(n):
                            A[i][j] = 4*p**2*(1.43e-34*(x[i][1] - x[j][1])**4 - \
                                    1.11e-27*(x[i][1] - x[j][1])**2 + 7.17e-22)*\
                                    np.exp(-3.87e-7*(x[i][0] - x[j][0])**2 - \
                                    3.87e-7*(x[i][1] - x[j][1])**2)
                elif mode == 'cross':
                    A = np.zeros((n,nn))
                    for i in range(n):
                        for j in range(nn):
                            A[i][j] = 4*p**2*(1.43e-34*(x[i][1] - z[j][1])**4 - \
                                    1.11e-27*(x[i][1] - z[j][1])**2 + 7.17e-22)*\
                                    np.exp(-3.87e-7*(x[i][0] - z[j][0])**2 - \
                                    3.87e-7*(x[i][1] - z[j][1])**2)
            else:
                raise Exception("Calling a derivative in RBF that doesn't exist")
            return A
except:
    pass
In [2]:
# Has to be added to pyGPs/Testing/unit_test_cov.py and then unit_test_cov.py
# has to be executed.
# For testing purposes only:

def test_cov_new(self):
    k = pyGPs.cov.MyKernel()     # specify your covariance function
    self.checkCovariance(k)