From 8ef7706fe557924c721c572b1f766530b7fbf10b Mon Sep 17 00:00:00 2001 From: Mark Pullin Date: Mon, 2 Jul 2018 15:23:58 +0100 Subject: [PATCH 1/4] Add symmetric kernel --- GPy/kern/__init__.py | 1 + GPy/kern/src/symmetric.py | 137 ++++++++++++++++++++++++++++++++++++ GPy/testing/kernel_tests.py | 14 +++- 3 files changed, 151 insertions(+), 1 deletion(-) create mode 100644 GPy/kern/src/symmetric.py diff --git a/GPy/kern/__init__.py b/GPy/kern/__init__.py index 96abab39..1fedb314 100644 --- a/GPy/kern/__init__.py +++ b/GPy/kern/__init__.py @@ -34,6 +34,7 @@ from .src.splitKern import DEtime as DiffGenomeKern from .src.spline import Spline from .src.basis_funcs import LogisticBasisFuncKernel, LinearSlopeBasisFuncKernel, BasisFuncKernel, ChangePointBasisFuncKernel, DomainKernel, PolynomialBasisFuncKernel from .src.grid_kerns import GridRBF +from .src.symmetric import Symmetric from .src.sde_matern import sde_Matern32 from .src.sde_matern import sde_Matern52 diff --git a/GPy/kern/src/symmetric.py b/GPy/kern/src/symmetric.py new file mode 100644 index 00000000..3072b38c --- /dev/null +++ b/GPy/kern/src/symmetric.py @@ -0,0 +1,137 @@ +import numpy as np + +from .kern import Kern + + +class Symmetric(Kern): + """ + Symmetric kernel that models a function with even or odd symmetry: + + .. math:: + + f(x) = f(Ax) + + or + + .. math:: + + f(x) = -f(Ax) + + it does this by modelling: + + .. math:: + + f(x) = g(x) \pm g(Ax) + + with kernel + + .. math:: + + k(x, x') \pm k(Ax, x') \pm k(x, Ax') + k(Ax, Ax') + + where k(x, x') is the kernel of g(x) + """ + + def __init__(self, base_kernel, transform, symmetry_type): + + super().__init__(1, [0], name='symmetric_kernel') + if symmetry_type is 'odd': + self.symmetry_sign = -1. + elif symmetry_type is 'even': + self.symmetry_sign = 1. + else: + raise ValueError('symmetry_type input must be ''odd'' or ''even''') + self.transform = transform + self.base_kernel = base_kernel + self.param_names = base_kernel.parameter_names() + self.link_parameters(self.base_kernel) + + def K(self, X, X2): + X_sym = X.dot(self.transform) + + if X2 is None: + X2 = X + X2_sym = X_sym + else: + X2_sym = X2.dot(self.transform) + + cross_term_x_ax = self.symmetry_sign * self.base_kernel.K(X, X2_sym) + + if X2 is None: + cross_term_ax_x = cross_term_x_ax.T + else: + cross_term_ax_x = self.symmetry_sign * \ + self.base_kernel.K(X_sym, X2) + + return (self.base_kernel.K(X, X2) + cross_term_x_ax + cross_term_ax_x + + self.base_kernel.K(X_sym, X2_sym)) + + def Kdiag(self, X): + n_points = X.shape[0] + X_sym = X.dot(self.transform) + + # Evaluate cross terms in batches, taking the diag of a larger matrix + # is wasteful, but is more efficient than calling kernel.K for each data point + batch_size = 100 + n_batches = int(np.ceil(n_points / batch_size)) + cross_term = np.zeros(X.shape[0]) + for i in range(n_batches): + i_start = i * batch_size + i_end = np.min([(i + 1) * batch_size, n_points]) + cross_term[i_start:i_end] = np.diag(self.base_kernel.K( + X_sym[i_start:i_end, :], X[i_start:i_end, :])) + + return self.base_kernel.Kdiag(X) + 2 * self.symmetry_sign * cross_term + self.base_kernel.Kdiag(X_sym) + + def update_gradients_full(self, dL_dK, X, X2): + X_sym = X.dot(self.transform) + if X2 is None: + X2 = X + X2_sym = X2.dot(self.transform) + + # Get gradients from base kernel one term at a time + self.base_kernel.update_gradients_full(dL_dK, X_sym, X2) + gradient = self.symmetry_sign * self.base_kernel.gradient.copy() + + self.base_kernel.update_gradients_full(dL_dK, X, X2_sym) + gradient += self.symmetry_sign * self.base_kernel.gradient.copy() + + self.base_kernel.update_gradients_full(dL_dK, X_sym, X2_sym) + gradient += self.base_kernel.gradient.copy() + + self.base_kernel.update_gradients_full(dL_dK, X, X2) + gradient += self.base_kernel.gradient.copy() + + # Set gradients + self.base_kernel.gradient = gradient + + def update_gradients_diag(self, dL_dK, X): + + dL_dK_full = np.diag(dL_dK) + X_sym = X.dot(self.transform) + + self.base_kernel.update_gradients_diag(dL_dK, X_sym) + gradient = self.base_kernel.gradient.copy() + + self.base_kernel.update_gradients_diag(dL_dK, X) + gradient += self.base_kernel.gradient.copy() + + # The contribution from both cross terms is the same + self.base_kernel.update_gradients_full(dL_dK_full, X, X_sym) + gradient += 2 * self.symmetry_sign * self.base_kernel.gradient.copy() + + self.base_kernel.gradient = gradient + + def gradients_X(self, dL_dK, X, X2): + X_sym = X.dot(self.transform) + if X2 is None: + X2 = X + X2_sym = X.dot(self.transform) + dL_dK = dL_dK + dL_dK.T + else: + X2_sym = X2.dot(self.transform) + + return (self.base_kernel.gradients_X(dL_dK, X, X2) + + self.base_kernel.gradients_X(dL_dK, X_sym, X2_sym).dot(self.transform.T) + + self.symmetry_sign * self.base_kernel.gradients_X(dL_dK, X, X2_sym) + + self.symmetry_sign * self.base_kernel.gradients_X(dL_dK, X_sym, X2).dot(self.transform.T)) diff --git a/GPy/testing/kernel_tests.py b/GPy/testing/kernel_tests.py index e1c9d934..02186c62 100644 --- a/GPy/testing/kernel_tests.py +++ b/GPy/testing/kernel_tests.py @@ -482,7 +482,19 @@ class KernelGradientTestsContinuous(unittest.TestCase): k = GPy.kern.StdPeriodic(self.D) k.randomize() self.assertTrue(check_kernel_gradient_functions(k, X=self.X, X2=self.X2, verbose=verbose)) - + + def test_symmetric_even(self): + k_base = GPy.kern.Linear(1) + GPy.kern.RBF(1) + transform = -np.array([[1.0]]) + k = GPy.kern.Symmetric(k_base, transform, 'even') + self.assertTrue(check_kernel_gradient_functions(k)) + + def test_symmetric_odd(self): + k_base = GPy.kern.Linear(1) + GPy.kern.RBF(1) + transform = -np.array([[1.0]]) + k = GPy.kern.Symmetric(k_base, transform, 'odd') + self.assertTrue(check_kernel_gradient_functions(k)) + def test_MultioutputKern(self): k1 = GPy.kern.RBF(self.D, ARD=True) k1.randomize() From 1a42cbd3426cebfac67935202621021c48b747b3 Mon Sep 17 00:00:00 2001 From: Mark Pullin Date: Mon, 2 Jul 2018 16:41:50 +0100 Subject: [PATCH 2/4] Add param descriptions --- GPy/kern/src/symmetric.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/GPy/kern/src/symmetric.py b/GPy/kern/src/symmetric.py index 3072b38c..1577d6bd 100644 --- a/GPy/kern/src/symmetric.py +++ b/GPy/kern/src/symmetric.py @@ -30,9 +30,13 @@ class Symmetric(Kern): k(x, x') \pm k(Ax, x') \pm k(x, Ax') + k(Ax, Ax') where k(x, x') is the kernel of g(x) + + :param base_kernel: kernel to make symmetric + :param transform: transformation matrix describing symmetry plane, A in equations above + :param symmetry_type: 'odd' or 'even' depending on the symmetry needed """ - def __init__(self, base_kernel, transform, symmetry_type): + def __init__(self, base_kernel, transform, symmetry_type='even'): super().__init__(1, [0], name='symmetric_kernel') if symmetry_type is 'odd': From 09c72eeec56e1cfd51c277050b7241fc757ef643 Mon Sep 17 00:00:00 2001 From: Mark Pullin Date: Mon, 2 Jul 2018 21:03:18 +0100 Subject: [PATCH 3/4] Make symmetric kernel work with python 2.7 --- GPy/kern/src/symmetric.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/GPy/kern/src/symmetric.py b/GPy/kern/src/symmetric.py index 1577d6bd..85c0bc61 100644 --- a/GPy/kern/src/symmetric.py +++ b/GPy/kern/src/symmetric.py @@ -38,7 +38,7 @@ class Symmetric(Kern): def __init__(self, base_kernel, transform, symmetry_type='even'): - super().__init__(1, [0], name='symmetric_kernel') + super(Symmetric, self).__init__(1, [0], name='symmetric_kernel') if symmetry_type is 'odd': self.symmetry_sign = -1. elif symmetry_type is 'even': @@ -77,7 +77,7 @@ class Symmetric(Kern): # Evaluate cross terms in batches, taking the diag of a larger matrix # is wasteful, but is more efficient than calling kernel.K for each data point batch_size = 100 - n_batches = int(np.ceil(n_points / batch_size)) + n_batches = int(np.ceil(n_points / float(batch_size))) cross_term = np.zeros(X.shape[0]) for i in range(n_batches): i_start = i * batch_size From 1cd860645131b52c86c7e61de78864213ebec7ae Mon Sep 17 00:00:00 2001 From: Mark Pullin Date: Thu, 19 Jul 2018 10:17:49 +0100 Subject: [PATCH 4/4] Expand class description and some speed improvements --- GPy/kern/src/symmetric.py | 45 ++++++++++++++++++++++++++++++++------- 1 file changed, 37 insertions(+), 8 deletions(-) diff --git a/GPy/kern/src/symmetric.py b/GPy/kern/src/symmetric.py index 85c0bc61..c7207023 100644 --- a/GPy/kern/src/symmetric.py +++ b/GPy/kern/src/symmetric.py @@ -7,11 +7,25 @@ class Symmetric(Kern): """ Symmetric kernel that models a function with even or odd symmetry: + For even symmetry we have: + .. math:: f(x) = f(Ax) - or + we then model the function as: + + .. math:: + + f(x) = g(x) + g(Ax) + + the corresponding kernel is: + + .. math:: + + k(x, x') + k(Ax, x') + k(x, Ax') + k(Ax, Ax') + + For odd symmetry we have: .. math:: @@ -21,13 +35,13 @@ class Symmetric(Kern): .. math:: - f(x) = g(x) \pm g(Ax) + f(x) = g(x) - g(Ax) with kernel .. math:: - k(x, x') \pm k(Ax, x') \pm k(x, Ax') + k(Ax, Ax') + k(x, x') - k(Ax, x') - k(x, Ax') + k(Ax, Ax') where k(x, x') is the kernel of g(x) @@ -37,8 +51,8 @@ class Symmetric(Kern): """ def __init__(self, base_kernel, transform, symmetry_type='even'): - - super(Symmetric, self).__init__(1, [0], name='symmetric_kernel') + n_dims = max(base_kernel.active_dims) + 1 + super(Symmetric, self).__init__(n_dims, list(range(n_dims)), name='symmetric_kernel') if symmetry_type is 'odd': self.symmetry_sign = -1. elif symmetry_type is 'even': @@ -114,15 +128,30 @@ class Symmetric(Kern): dL_dK_full = np.diag(dL_dK) X_sym = X.dot(self.transform) + # Calculate gradient for k(Ax, Ax') self.base_kernel.update_gradients_diag(dL_dK, X_sym) gradient = self.base_kernel.gradient.copy() + # Calculate gradient for k(x, x') self.base_kernel.update_gradients_diag(dL_dK, X) gradient += self.base_kernel.gradient.copy() - # The contribution from both cross terms is the same - self.base_kernel.update_gradients_full(dL_dK_full, X, X_sym) - gradient += 2 * self.symmetry_sign * self.base_kernel.gradient.copy() + # Batch process cross term for speed + batch_size = 100 + n_points = dL_dK.shape[0] + n_batches = int(np.ceil(n_points / float(batch_size))) + gradient_part = np.zeros(gradient.shape) + for i in range(n_batches): + i_start = i * batch_size + i_end = np.min([(i + 1) * batch_size, n_points]) + dL_dK_part = dL_dK_full[i_start:i_end, i_start:i_end] + X_part = X[i_start:i_end, :] + X_sym_part = X_sym[i_start:i_end, :] + self.base_kernel.update_gradients_full( + dL_dK_part, X_part, X_sym_part) + gradient_part += self.base_kernel.gradient.copy() + + gradient += 2 * self.symmetry_sign * gradient_part self.base_kernel.gradient = gradient