diff --git a/GPy/kern/src/kern.py b/GPy/kern/src/kern.py index bac85d58..6a7aea19 100644 --- a/GPy/kern/src/kern.py +++ b/GPy/kern/src/kern.py @@ -185,6 +185,9 @@ class Kern(Parameterized): def update_gradients_full(self, dL_dK, X, X2): """Set the gradients of all parameters when doing full (N) inference.""" raise NotImplementedError + + def reset_gradients(self): + raise NotImplementedError def update_gradients_expectations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior): """ @@ -348,7 +351,7 @@ class CombinationKernel(Kern): A combination kernel combines (a list of) kernels and works on those. Examples are the HierarchicalKernel or Add and Prod kernels. """ - def __init__(self, kernels, name, extra_dims=[]): + def __init__(self, kernels, name, extra_dims=[], link_parameters=True): """ Abstract super class for combination kernels. A combination kernel combines (a list of) kernels and works on those. @@ -372,7 +375,8 @@ class CombinationKernel(Kern): self._all_dims_active = np.array(np.concatenate((np.arange(effective_input_dim), extra_dims if extra_dims is not None else [])), dtype=int) self.extra_dims = extra_dims - self.link_parameters(*kernels) + if link_parameters: + self.link_parameters(*kernels) def _to_dict(self): input_dict = super(CombinationKernel, self)._to_dict() diff --git a/GPy/kern/src/kernel_slice_operations.py b/GPy/kern/src/kernel_slice_operations.py index 327436e7..50a1e2e8 100644 --- a/GPy/kern/src/kernel_slice_operations.py +++ b/GPy/kern/src/kernel_slice_operations.py @@ -97,17 +97,17 @@ def _slice_Kdiag(f): def _slice_update_gradients_full(f): @wraps(f) - def wrap(self, dL_dK, X, X2=None): + def wrap(self, dL_dK, X, X2=None, *a, **kw): with _Slice_wrap(self, X, X2) as s: - ret = f(self, dL_dK, s.X, s.X2) + ret = f(self, dL_dK, s.X, s.X2, *a, **kw) return ret return wrap def _slice_update_gradients_diag(f): @wraps(f) - def wrap(self, dL_dKdiag, X): + def wrap(self, dL_dKdiag, X, *a, **kw): with _Slice_wrap(self, X, None) as s: - ret = f(self, dL_dKdiag, s.X) + ret = f(self, dL_dKdiag, s.X, *a, **kw) return ret return wrap diff --git a/GPy/kern/src/multioutput_kern.py b/GPy/kern/src/multioutput_kern.py index 4614ef5a..c196e2cd 100644 --- a/GPy/kern/src/multioutput_kern.py +++ b/GPy/kern/src/multioutput_kern.py @@ -1,11 +1,11 @@ from .kern import Kern, CombinationKernel import numpy as np from functools import reduce, partial -from GPy.util.multioutput import index_to_slices +from .independent_outputs import index_to_slices from paramz.caching import Cache_this class MultioutputKern(CombinationKernel): - def __init__(self, kernels, cross_covariances, name='MultioutputKern'): + def __init__(self, kernels, cross_covariances={}, name='MultioutputKern'): #kernels contains a list of kernels as input, if not isinstance(kernels, list): self.single_kern = True @@ -20,7 +20,7 @@ class MultioutputKern(CombinationKernel): # after slicing. This is why the index_dim is just the last column: self.index_dim = -1 - super(MultioutputKern, self).__init__(kernels=kernels, extra_dims=[self.index_dim], name=name, link_params=False) + super(MultioutputKern, self).__init__(kernels=kernels, extra_dims=[self.index_dim], name=name, link_parameters=False) nl = len(kernels) #build covariance structure @@ -36,7 +36,7 @@ class MultioutputKern(CombinationKernel): elif cross_covariances.get((i,j)) is not None: #cross covariance is given covariance[i][j] = cross_covariances.get((i,j)) else: # zero matrix - covariance[i][j] = {'K': lambda x, x2: np.zeros((x.shape[0],x2.shape[0])), 'update_gradients_full': lambda x, x2: np.zeros((x.shape[0],x2.shape[0])), 'gradients_X': lambda x, x2: np.zeros((x.shape[0],x.shape[1]))} + covariance[i][j] = {'K': lambda x, x2: np.zeros((x.shape[0],x2.shape[0])), 'update_gradients_full': lambda dl_dk, x, x2, reset: np.zeros(dl_dk.shape), 'gradients_X': lambda dl_dk, x, x2: np.zeros((x.shape[0],x.shape[1]))} if unique is True: linked.append(i) self.covariance = covariance @@ -63,26 +63,27 @@ class MultioutputKern(CombinationKernel): def reset_gradients(self): for kern in self.kern: kern.reset_gradients() - def update_gradients_full(self,dL_dK,X,X2=None, reset=True): - if reset: - self.reset_gradients() - if X2 is None: - X2 = X + def update_gradients_full(self,dL_dK, X, X2=None, reset=True): + if reset: self.reset_gradients() slices = index_to_slices(X[:,self.index_dim]) - slices2 = index_to_slices(X2[:,self.index_dim]) - [[[[ self.covariance[i][j]['update_gradients_full'](dL_dK[slices[i][k],slices2[j][l]], X[slices[i][k],:], X2[slices2[j][l],:], False) for k in range(len(slices[i]))] for l in range(len(slices2[j]))] for i in range(len(slices))] for j in range(len(slices2))] - - def update_gradients_diag(self, dL_dKdiag, X): - for kern in self.kerns: kern.reset_gradients() + if X2 is not None: + slices2 = index_to_slices(X2[:,self.index_dim]) + [[[[ self.covariance[i][j]['update_gradients_full'](dL_dK[slices[i][k],slices2[j][l]], X[slices[i][k],:], X2[slices2[j][l],:], False) for k in range(len(slices[i]))] for l in range(len(slices2[j]))] for i in range(len(slices))] for j in range(len(slices2))] + else: + [[[[ self.covariance[i][j]['update_gradients_full'](dL_dK[slices[i][k],slices[j][l]], X[slices[i][k],:], X[slices[j][l],:] , False) for k in range(len(slices[i]))] for l in range(len(slices[j]))] for i in range(len(slices))] for j in range(len(slices))] + + def update_gradients_diag(self, dL_dKdiag, X, reset=True): + if reset: self.reset_gradients() slices = index_to_slices(X[:,self.index_dim]) kerns = itertools.repeat(self.kern) if self.single_kern else self.kern - [[ self.kerns[i].update_gradients_diag(dL_dKdiag[slices[i][k]], X[slices[i][k],:], False) for k in range(len(slices[i]))] for i in range(len(slices))] + [[ self.kern[i].update_gradients_diag(dL_dKdiag[slices[i][k]], X[slices[i][k],:], False) for k in range(len(slices[i]))] for i in range(len(slices))] def gradients_X(self,dL_dK, X, X2=None): - if X2 is None: - X2 = X slices = index_to_slices(X[:,self.index_dim]) - slices2 = index_to_slices(X2[:,self.index_dim]) target = np.zeros((X.shape[0], X.shape[1]) ) - [[[[ target.__setitem__((slices[i][k]), target[slices[i][k],:] + self.covariance[i][j]['gradients_X'](dL_dK[slices[i][k],slices2[j][l]], X[slices[i][k],:], X2[slices2[j][l],:]) ) for k in range(len(slices[i]))] for l in range(len(slices2[j]))] for i in range(len(slices))] for j in range(len(slices2))] + if X2 is not None: + slices2 = index_to_slices(X2[:,self.index_dim]) + [[[[ target.__setitem__((slices[i][k]), target[slices[i][k],:] + self.covariance[i][j]['gradients_X'](dL_dK[slices[i][k],slices2[j][l]], X[slices[i][k],:], X2[slices2[j][l],:]) ) for k in range(len(slices[i]))] for l in range(len(slices2[j]))] for i in range(len(slices))] for j in range(len(slices2))] + else: + [[[[ target.__setitem__((slices[i][k]), target[slices[i][k],:] + self.covariance[i][j]['gradients_X'](dL_dK[slices[i][k],slices[j][l]], X[slices[i][k],:], (None if (i==j and k==l) else X[slices[j][l],:] )) ) for k in range(len(slices[i]))] for l in range(len(slices[j]))] for i in range(len(slices))] for j in range(len(slices))] return target \ No newline at end of file diff --git a/GPy/kern/src/rbf.py b/GPy/kern/src/rbf.py index 0b6730d8..479e1357 100644 --- a/GPy/kern/src/rbf.py +++ b/GPy/kern/src/rbf.py @@ -107,10 +107,10 @@ class RBF(Stationary): def gradients_qX_expectations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior): return self.psicomp.psiDerivativecomputations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior)[3:] - def update_gradients_diag(self, dL_dKdiag, X): - super(RBF,self).update_gradients_diag(dL_dKdiag, X) + def update_gradients_diag(self, dL_dKdiag, X, reset=True): + super(RBF,self).update_gradients_diag(dL_dKdiag, X, reset=reset) if self.use_invLengthscale: self.inv_l.gradient =self.lengthscale.gradient*(self.lengthscale**3/-2.) - def update_gradients_full(self, dL_dK, X, X2=None): - super(RBF,self).update_gradients_full(dL_dK, X, X2) + def update_gradients_full(self, dL_dK, X, X2=None, reset=True): + super(RBF,self).update_gradients_full(dL_dK, X, X2, reset=reset) if self.use_invLengthscale: self.inv_l.gradient =self.lengthscale.gradient*(self.lengthscale**3/-2.) diff --git a/GPy/kern/src/stationary.py b/GPy/kern/src/stationary.py index 4e8ddb77..cd09a4a2 100644 --- a/GPy/kern/src/stationary.py +++ b/GPy/kern/src/stationary.py @@ -171,7 +171,14 @@ class Stationary(Kern): ret[:] = self.variance return ret - def update_gradients_diag(self, dL_dKdiag, X): + def reset_gradients(self): + self.variance.gradient = 0. + if not self.ARD: + self.lengthscale.gradient = 0. + else: + self.lengthscale.gradient = np.zeros(self.input_dim) + + def update_gradients_diag(self, dL_dKdiag, X, reset=True): """ Given the derivative of the objective with respect to the diagonal of the covariance matrix, compute the derivative wrt the parameters of @@ -179,16 +186,18 @@ class Stationary(Kern): See also update_gradients_full """ - self.variance.gradient = np.sum(dL_dKdiag) - self.lengthscale.gradient = 0. + if reset: self.reset_gradients() + self.variance.gradient += np.sum(dL_dKdiag) + self.lengthscale.gradient += 0. - def update_gradients_full(self, dL_dK, X, X2=None): + def update_gradients_full(self, dL_dK, X, X2=None, reset=True): """ Given the derivative of the objective wrt the covariance matrix (dL_dK), compute the gradient wrt the parameters of this kernel, and store in the parameters object as e.g. self.variance.gradient """ - self.variance.gradient = np.sum(self.K(X, X2)* dL_dK)/self.variance + if reset: self.reset_gradients() + self.variance.gradient += np.sum(self.K(X, X2)* dL_dK)/self.variance #now the lengthscale gradient(s) dL_dr = self.dK_dr_via_X(X, X2) * dL_dK @@ -197,12 +206,12 @@ class Stationary(Kern): tmp = dL_dr*self._inv_dist(X, X2) if X2 is None: X2 = X if config.getboolean('cython', 'working'): - self.lengthscale.gradient = self._lengthscale_grads_cython(tmp, X, X2) + self.lengthscale.gradient += self._lengthscale_grads_cython(tmp, X, X2) else: - self.lengthscale.gradient = self._lengthscale_grads_pure(tmp, X, X2) + self.lengthscale.gradient += self._lengthscale_grads_pure(tmp, X, X2) else: r = self._scaled_dist(X, X2) - self.lengthscale.gradient = -np.sum(dL_dr*r)/self.lengthscale + self.lengthscale.gradient += -np.sum(dL_dr*r)/self.lengthscale def update_gradients_direct(self, dL_dVar, dL_dLen): @@ -632,7 +641,7 @@ class RatQuad(Stationary): def dK_dr(self, r): r2 = np.square(r) # return -self.variance*self.power*r*np.power(1. + r2/2., - self.power - 1.) - return-self.variance*self.power*r*np.exp(-(self.power+1)*np.log1p(r2/2.)) + return -self.variance*self.power*r*np.exp(-(self.power+1)*np.log1p(r2/2.)) def update_gradients_full(self, dL_dK, X, X2=None): super(RatQuad, self).update_gradients_full(dL_dK, X, X2) diff --git a/GPy/testing/kernel_tests.py b/GPy/testing/kernel_tests.py index 0d4db63b..f903c9aa 100644 --- a/GPy/testing/kernel_tests.py +++ b/GPy/testing/kernel_tests.py @@ -484,11 +484,16 @@ class KernelGradientTestsContinuous(unittest.TestCase): self.assertTrue(check_kernel_gradient_functions(k, X=self.X, X2=self.X2, verbose=verbose)) def test_MultioutputKern(self): - k1 = GPy.kern.RBF(self.D-1, ARD=True) + k1 = GPy.kern.RBF(self.D, ARD=True) k1.randomize() - k2 = GPy.kern.RBF(self.D-1, ARD=True) + k2 = GPy.kern.RBF(self.D, ARD=True) k2.randomize() - k = GPy.kern.MultioutputKern([k1,k2],) + + k = GPy.kern.MultioutputKern([k1, k2]) + Xt,_,_ = GPy.util.multioutput.build_XY([self.X, self.X]) + X2t,_,_ = GPy.util.multioutput.build_XY([self.X2, self.X2]) + self.assertTrue(check_kernel_gradient_functions(k, X=Xt, X2=X2t, verbose=verbose, fixed_X_dims=-1)) + def test_Precomputed(self): Xall = np.concatenate([self.X, self.X2])