From a9c8ef817af5e85b656fb33d83540fbf42d2188f Mon Sep 17 00:00:00 2001 From: alessandratosi Date: Wed, 20 Apr 2016 12:09:25 +0100 Subject: [PATCH 01/20] Update function kern.gradients_XX() to compute cross-covariance terms --- GPy/kern/src/add.py | 21 ++++++--- GPy/kern/src/kern.py | 2 +- GPy/kern/src/kernel_slice_operations.py | 31 ++++++++------ GPy/kern/src/stationary.py | 57 ++++++++++++++++--------- 4 files changed, 69 insertions(+), 42 deletions(-) diff --git a/GPy/kern/src/add.py b/GPy/kern/src/add.py index ec09f157..2ad515dd 100644 --- a/GPy/kern/src/add.py +++ b/GPy/kern/src/add.py @@ -85,13 +85,20 @@ class Add(CombinationKernel): [target.__iadd__(p.gradients_X_diag(dL_dKdiag, X)) for p in self.parts] return target - def gradients_XX(self, dL_dK, X, X2): - if X2 is None: - target = np.zeros((X.shape[0], X.shape[0], X.shape[1])) - else: - target = np.zeros((X.shape[0], X2.shape[0], X.shape[1])) - [target.__iadd__(p.gradients_XX(dL_dK, X, X2)) for p in self.parts] - return target + # def gradients_XX(self, dL_dK, X, X2, cov=True): + # if cov==True: # full covarance + # if X2 is None: + # target = np.zeros((X.shape[0], X.shape[0], X.shape[1], X.shape[1])) + # else: + # target = np.zeros((X.shape[0], X2.shape[0], X.shape[1], X.shape[1])) + # else: # diagonal covariance + # if X2 is None: + # target = np.zeros((X.shape[0], X.shape[0], X.shape[1])) + # else: + # target = np.zeros((X.shape[0], X2.shape[0], X.shape[1])) + + # [target.__iadd__(p.gradients_XX(dL_dK, X, X2, cov=True)) for p in self.parts] + # return target def gradients_XX_diag(self, dL_dKdiag, X): target = np.zeros(X.shape) diff --git a/GPy/kern/src/kern.py b/GPy/kern/src/kern.py index 393cf1c2..b64d145b 100644 --- a/GPy/kern/src/kern.py +++ b/GPy/kern/src/kern.py @@ -132,7 +132,7 @@ class Kern(Parameterized): raise NotImplementedError def gradients_X_X2(self, dL_dK, X, X2): return self.gradients_X(dL_dK, X, X2), self.gradients_X(dL_dK.T, X2, X) - def gradients_XX(self, dL_dK, X, X2): + def gradients_XX(self, dL_dK, X, X2, cov='False'): """ .. math:: diff --git a/GPy/kern/src/kernel_slice_operations.py b/GPy/kern/src/kernel_slice_operations.py index 57b34de9..c0d46c0f 100644 --- a/GPy/kern/src/kernel_slice_operations.py +++ b/GPy/kern/src/kernel_slice_operations.py @@ -24,7 +24,7 @@ class KernCallsViaSlicerMeta(ParametersChangedMeta): put_clean(dct, 'update_gradients_diag', _slice_update_gradients_diag) put_clean(dct, 'gradients_X', _slice_gradients_X) put_clean(dct, 'gradients_X_X2', _slice_gradients_X) - put_clean(dct, 'gradients_XX', _slice_gradients_XX) +# put_clean(dct, 'gradients_XX', _slice_gradients_XX) put_clean(dct, 'gradients_XX_diag', _slice_gradients_X_diag) put_clean(dct, 'gradients_X_diag', _slice_gradients_X_diag) @@ -112,18 +112,23 @@ def _slice_gradients_X(f): return ret return wrap -def _slice_gradients_XX(f): - @wraps(f) - def wrap(self, dL_dK, X, X2=None): - if X2 is None: - N, M = X.shape[0], X.shape[0] - else: - N, M = X.shape[0], X2.shape[0] - with _Slice_wrap(self, X, X2, ret_shape=(N, M, X.shape[1])) as s: - #with _Slice_wrap(self, X, X2, ret_shape=None) as s: - ret = s.handle_return_array(f(self, dL_dK, s.X, s.X2)) - return ret - return wrap +# def _slice_gradients_XX(f): +# @wraps(f) +# def wrap(self, dL_dK, X, X2=None, cov=True): +# if X2 is None: +# N, M = X.shape[0], X.shape[0] +# else: +# N, M = X.shape[0], X2.shape[0] +# if cov==True: # full covariance +# with _Slice_wrap(self, X, X2, ret_shape=(N, M, X.shape[1], X.shape[1])) as s: +# #with _Slice_wrap(self, X, X2, ret_shape=None) as s: +# ret = s.handle_return_array(f(self, dL_dK, s.X, s.X2)) +# else: # diagonal covariance +# with _Slice_wrap(self, X, X2, ret_shape=(N, M, X.shape[1])) as s: +# #with _Slice_wrap(self, X, X2, ret_shape=None) as s: +# ret = s.handle_return_array(f(self, dL_dK, s.X, s.X2)) +# return ret +# return wrap def _slice_gradients_X_diag(f): @wraps(f) diff --git a/GPy/kern/src/stationary.py b/GPy/kern/src/stationary.py index 5e137abb..90aa4297 100644 --- a/GPy/kern/src/stationary.py +++ b/GPy/kern/src/stationary.py @@ -218,45 +218,60 @@ class Stationary(Kern): else: return self._gradients_X_pure(dL_dK, X, X2) - def gradients_XX(self, dL_dK, X, X2=None): + def gradients_XX(self, dL_dK, X, X2=None, cov=True): """ Given the derivative of the objective K(dL_dK), compute the second derivative of K wrt X and X2: + cov = Full: returns the full covariance matrix [QxQ] of the input dimensionfor each pair or vectors + cov = Diag: returns the diagonal of the covariance matrix [QxQ] of the input dimensionfor each pair + or vectors (computationally more efficient if the full covariance matrix is not needed) ..math: - \frac{\partial^2 K}{\partial X\partial X2} + \frac{\partial^2 K}{\partial X2 ^2} = - \frac{\partial^2 K}{\partial X\partial X2} ..returns: - dL2_dXdX2: NxMxQ, for X [NxQ] and X2[MxQ] (X2 is X if, X2 is None) - Thus, we return the second derivative in X2. + dL2_dXdX2: [NxMxQ] in the cov=Diag case, or [NxMxQxQ] in the cov=full case, + for X [NxQ] and X2[MxQ] (X2 is X if, X2 is None) + Thus, we return the second derivative in X2. """ - # The off diagonals in Q are always zero, this should also be true for the Linear kernel... # According to multivariable chain rule, we can chain the second derivative through r: # d2K_dXdX2 = dK_dr*d2r_dXdX2 + d2K_drdr * dr_dX * dr_dX2: invdist = self._inv_dist(X, X2) invdist2 = invdist**2 - - dL_dr = self.dK_dr_via_X(X, X2) * dL_dK + dL_dr = self.dK_dr_via_X(X, X2) # * dL_dK we perofrm this product later tmp1 = dL_dr * invdist - - dL_drdr = self.dK2_drdr_via_X(X, X2) * dL_dK + dL_drdr = self.dK2_drdr_via_X(X, X2) # * dL_dK we perofrm this product later tmp2 = dL_drdr * invdist2 - - l2 = np.ones(X.shape[1]) * self.lengthscale**2 + l2 = np.ones(X.shape[1])*self.lengthscale**2 #np.multiply(np.ones(X.shape[1]) ,self.lengthscale**2) + print ['l2',l2] if X2 is None: X2 = X tmp1 -= np.eye(X.shape[0])*self.variance else: - tmp1[X==X2.T] -= self.variance - - grad = np.empty((X.shape[0], X2.shape[0], X.shape[1]), dtype=np.float64) - #grad = np.empty(X.shape, dtype=np.float64) - for q in range(self.input_dim): - tmpdist2 = (X[:,[q]]-X2[:,[q]].T) ** 2 - grad[:, :, q] = ((tmp1*invdist2 - tmp2)*tmpdist2/l2[q] - tmp1)/l2[q] - #grad[:, :, q] = ((tmp1*(((tmpdist2)*invdist2/l2[q])-1)) - (tmp2*(tmpdist2))/l2[q])/l2[q] - #np.sum(((tmp1*(((tmpdist2)*invdist2/l2[q])-1)) - (tmp2*(tmpdist2))/l2[q])/l2[q], axis=1, out=grad[:,q]) - #np.sum( - (tmp2*(tmpdist**2)), axis=1, out=grad[:,q]) + #tmp1[X==X2.T] -= self.variance # Old version, to be removed + # (seems to have a bug: it is subtracted to the first X1 anyway) + tmp1[invdist2==0.] -= self.variance + + if cov==True: # full covariance + grad = np.empty((X.shape[0], X2.shape[0], X2.shape[1], X.shape[1]), dtype=np.float64) + for q in range(self.input_dim): + for r in range(self.input_dim): + tmpdist2 = (X[:,[q]]-X2[:,[q]].T)*(X[:,[r]]-X2[:,[r]].T) # Introduce temporary distance + if r==q: + grad[:, :, q, r] = np.multiply(dL_dK,(np.multiply((tmp1*invdist2 - tmp2),tmpdist2)/l2[r] - tmp1)/l2[q]) + else: + grad[:, :, q, r] = np.multiply(dL_dK,(np.multiply((tmp1*invdist2 - tmp2),tmpdist2)/l2[r])/l2[q]) + else: + # Diagonal covariance + grad = np.empty((X.shape[0], X2.shape[0], X.shape[1]), dtype=np.float64) + #grad = np.empty(X.shape, dtype=np.float64) + for q in range(self.input_dim): + tmpdist2 = (X[:,[q]]-X2[:,[q]].T) ** 2 + grad[:, :, q] = np.multiply(dL_dK,(np.multiply((tmp1*invdist2 - tmp2),tmpdist2)/l2[q] - tmp1)/l2[q]) + #grad[:, :, q] = ((tmp1*invdist2 - tmp2)*tmpdist2/l2[q] - tmp1)/l2[q] + #grad[:, :, q] = ((tmp1*(((tmpdist2)*invdist2/l2[q])-1)) - (tmp2*(tmpdist2))/l2[q])/l2[q] + #np.sum(((tmp1*(((tmpdist2)*invdist2/l2[q])-1)) - (tmp2*(tmpdist2))/l2[q])/l2[q], axis=1, out=grad[:,q]) + #np.sum( - (tmp2*(tmpdist**2)), axis=1, out=grad[:,q]) return grad def gradients_XX_diag(self, dL_dK, X): From a1e4728f8a63774dc372be2f2651a7049b25d47d Mon Sep 17 00:00:00 2001 From: alessandratosi Date: Wed, 20 Apr 2016 14:50:39 +0100 Subject: [PATCH 02/20] added kernel tests for gradients_XX --- GPy/kern/src/add.py | 26 +++++------ GPy/kern/src/kern.py | 2 +- GPy/kern/src/kernel_slice_operations.py | 36 ++++++++-------- GPy/kern/src/stationary.py | 3 +- GPy/testing/kernel_tests.py | 57 +++++++++++++++++++++++++ 5 files changed, 90 insertions(+), 34 deletions(-) diff --git a/GPy/kern/src/add.py b/GPy/kern/src/add.py index 2ad515dd..9b83d633 100644 --- a/GPy/kern/src/add.py +++ b/GPy/kern/src/add.py @@ -85,20 +85,20 @@ class Add(CombinationKernel): [target.__iadd__(p.gradients_X_diag(dL_dKdiag, X)) for p in self.parts] return target - # def gradients_XX(self, dL_dK, X, X2, cov=True): - # if cov==True: # full covarance - # if X2 is None: - # target = np.zeros((X.shape[0], X.shape[0], X.shape[1], X.shape[1])) - # else: - # target = np.zeros((X.shape[0], X2.shape[0], X.shape[1], X.shape[1])) - # else: # diagonal covariance - # if X2 is None: - # target = np.zeros((X.shape[0], X.shape[0], X.shape[1])) - # else: - # target = np.zeros((X.shape[0], X2.shape[0], X.shape[1])) + def gradients_XX(self, dL_dK, X, X2, cov=True): + if cov==True: # full covarance + if X2 is None: + target = np.zeros((X.shape[0], X.shape[0], X.shape[1], X.shape[1])) + else: + target = np.zeros((X.shape[0], X2.shape[0], X.shape[1], X.shape[1])) + else: # diagonal covariance + if X2 is None: + target = np.zeros((X.shape[0], X.shape[0], X.shape[1])) + else: + target = np.zeros((X.shape[0], X2.shape[0], X.shape[1])) - # [target.__iadd__(p.gradients_XX(dL_dK, X, X2, cov=True)) for p in self.parts] - # return target + [target.__iadd__(p.gradients_XX(dL_dK, X, X2, cov)) for p in self.parts] + return target def gradients_XX_diag(self, dL_dKdiag, X): target = np.zeros(X.shape) diff --git a/GPy/kern/src/kern.py b/GPy/kern/src/kern.py index b64d145b..37307e6b 100644 --- a/GPy/kern/src/kern.py +++ b/GPy/kern/src/kern.py @@ -132,7 +132,7 @@ class Kern(Parameterized): raise NotImplementedError def gradients_X_X2(self, dL_dK, X, X2): return self.gradients_X(dL_dK, X, X2), self.gradients_X(dL_dK.T, X2, X) - def gradients_XX(self, dL_dK, X, X2, cov='False'): + def gradients_XX(self, dL_dK, X, X2, cov='True'): """ .. math:: diff --git a/GPy/kern/src/kernel_slice_operations.py b/GPy/kern/src/kernel_slice_operations.py index c0d46c0f..ddb16ea1 100644 --- a/GPy/kern/src/kernel_slice_operations.py +++ b/GPy/kern/src/kernel_slice_operations.py @@ -24,7 +24,7 @@ class KernCallsViaSlicerMeta(ParametersChangedMeta): put_clean(dct, 'update_gradients_diag', _slice_update_gradients_diag) put_clean(dct, 'gradients_X', _slice_gradients_X) put_clean(dct, 'gradients_X_X2', _slice_gradients_X) -# put_clean(dct, 'gradients_XX', _slice_gradients_XX) + put_clean(dct, 'gradients_XX', _slice_gradients_XX) put_clean(dct, 'gradients_XX_diag', _slice_gradients_X_diag) put_clean(dct, 'gradients_X_diag', _slice_gradients_X_diag) @@ -112,23 +112,23 @@ def _slice_gradients_X(f): return ret return wrap -# def _slice_gradients_XX(f): -# @wraps(f) -# def wrap(self, dL_dK, X, X2=None, cov=True): -# if X2 is None: -# N, M = X.shape[0], X.shape[0] -# else: -# N, M = X.shape[0], X2.shape[0] -# if cov==True: # full covariance -# with _Slice_wrap(self, X, X2, ret_shape=(N, M, X.shape[1], X.shape[1])) as s: -# #with _Slice_wrap(self, X, X2, ret_shape=None) as s: -# ret = s.handle_return_array(f(self, dL_dK, s.X, s.X2)) -# else: # diagonal covariance -# with _Slice_wrap(self, X, X2, ret_shape=(N, M, X.shape[1])) as s: -# #with _Slice_wrap(self, X, X2, ret_shape=None) as s: -# ret = s.handle_return_array(f(self, dL_dK, s.X, s.X2)) -# return ret -# return wrap +def _slice_gradients_XX(f): + @wraps(f) + def wrap(self, dL_dK, X, X2=None, cov=True): + if X2 is None: + N, M = X.shape[0], X.shape[0] + else: + N, M = X.shape[0], X2.shape[0] + if cov==True: # full covariance + with _Slice_wrap(self, X, X2, ret_shape=(N, M, X.shape[1], X.shape[1])) as s: + #with _Slice_wrap(self, X, X2, ret_shape=None) as s: + ret = s.handle_return_array(f(self, dL_dK, s.X, s.X2, cov=True)) + else: # diagonal covariance + with _Slice_wrap(self, X, X2, ret_shape=(N, M, X.shape[1])) as s: + #with _Slice_wrap(self, X, X2, ret_shape=None) as s: + ret = s.handle_return_array(f(self, dL_dK, s.X, s.X2, cov=True)) + return ret + return wrap def _slice_gradients_X_diag(f): @wraps(f) diff --git a/GPy/kern/src/stationary.py b/GPy/kern/src/stationary.py index 90aa4297..4de52d91 100644 --- a/GPy/kern/src/stationary.py +++ b/GPy/kern/src/stationary.py @@ -242,7 +242,6 @@ class Stationary(Kern): dL_drdr = self.dK2_drdr_via_X(X, X2) # * dL_dK we perofrm this product later tmp2 = dL_drdr * invdist2 l2 = np.ones(X.shape[1])*self.lengthscale**2 #np.multiply(np.ones(X.shape[1]) ,self.lengthscale**2) - print ['l2',l2] if X2 is None: X2 = X @@ -262,7 +261,7 @@ class Stationary(Kern): else: grad[:, :, q, r] = np.multiply(dL_dK,(np.multiply((tmp1*invdist2 - tmp2),tmpdist2)/l2[r])/l2[q]) else: - # Diagonal covariance + # Diagonal covariance, old code grad = np.empty((X.shape[0], X2.shape[0], X.shape[1]), dtype=np.float64) #grad = np.empty(X.shape, dtype=np.float64) for q in range(self.input_dim): diff --git a/GPy/testing/kernel_tests.py b/GPy/testing/kernel_tests.py index ae9aebfb..f47e9805 100644 --- a/GPy/testing/kernel_tests.py +++ b/GPy/testing/kernel_tests.py @@ -101,7 +101,21 @@ class Kern_check_dKdiag_dX(Kern_check_dK_dX): def parameters_changed(self): self.X.gradient[:] = self.kernel.gradients_X_diag(self.dL_dK.diagonal(), self.X) +class Kern_check_d2K_dXdX(Kern_check_model): + """This class allows gradient checks for the secondderivative of a kernel with respect to X. """ + def __init__(self, kernel=None, dL_dK=None, X=None, X2=None): + Kern_check_model.__init__(self,kernel=kernel,dL_dK=dL_dK, X=X, X2=X2) + self.X = Param('X',X) + self.link_parameter(self.X) + def log_likelihood(self): + return np.sum(self.kernel.gradients_X(self.dL_dK,self.X, self.X2)) + + def parameters_changed(self): + self.X.gradient[:] = self.kernel.gradients_XX(self.dL_dK, self.X, self.X2,cov=True).sum(0).sum(1) + +# class Kern_check_d2Kdiag_dXdX(Kern_check_model): +# """This class allows gradient checks for the secondderivative of a kernel diagonal with respect to X. """ def check_kernel_gradient_functions(kern, X=None, X2=None, output_ind=None, verbose=False, fixed_X_dims=None): """ @@ -239,6 +253,49 @@ def check_kernel_gradient_functions(kern, X=None, X2=None, output_ind=None, verb assert(result) return False + if verbose: + print("Checking gradients of dK(X, X) wrt X.") + try: + testmodel = Kern_check_d2K_dXdX(kern, X=X, X2=None) + if fixed_X_dims is not None: + testmodel.X[:,fixed_X_dims].fix() + result = testmodel.checkgrad(verbose=verbose) + except NotImplementedError: + result=True + if verbose: + print(("gradients_XX not implemented for " + kern.name)) + if result and verbose: + print("Check passed.") + if not result: + print(("Gradient of dK(X, X) wrt X failed for " + kern.name + " covariance function. Gradient values as follows:")) + testmodel.checkgrad(verbose=True) + assert(result) + pass_checks = False + return False + + if verbose: + print("Checking gradients of dK(X, X2) wrt X.") + try: + testmodel = Kern_check_d2K_dXdX(kern, X=X, X2=X2) + if fixed_X_dims is not None: + testmodel.X[:,fixed_X_dims].fix() + result = testmodel.checkgrad(verbose=verbose) + except NotImplementedError: + result=True + if verbose: + print(("gradients_XX not implemented for " + kern.name)) + if result and verbose: + print("Check passed.") + if not result: + print(("Gradient of dK(X, X2) wrt X failed for " + kern.name + " covariance function. Gradient values as follows:")) + testmodel.checkgrad(verbose=True) + assert(result) + pass_checks = False + return False + +# if verbose: +# print("Checking gradients of dKdiag(X, X) wrt X.") + return pass_checks From bbd8264235121780af99794a5a053ea7697ffc06 Mon Sep 17 00:00:00 2001 From: alessandratosi Date: Thu, 21 Apr 2016 12:09:39 +0100 Subject: [PATCH 03/20] modified kernel tests for gradients_XX --- GPy/testing/kernel_tests.py | 13 +++++-------- 1 file changed, 5 insertions(+), 8 deletions(-) diff --git a/GPy/testing/kernel_tests.py b/GPy/testing/kernel_tests.py index f47e9805..402599d2 100644 --- a/GPy/testing/kernel_tests.py +++ b/GPy/testing/kernel_tests.py @@ -105,14 +105,14 @@ class Kern_check_d2K_dXdX(Kern_check_model): """This class allows gradient checks for the secondderivative of a kernel with respect to X. """ def __init__(self, kernel=None, dL_dK=None, X=None, X2=None): Kern_check_model.__init__(self,kernel=kernel,dL_dK=dL_dK, X=X, X2=X2) - self.X = Param('X',X) - self.link_parameter(self.X) + #self.X = Param('X',X) + #self.link_parameter(self.X) def log_likelihood(self): return np.sum(self.kernel.gradients_X(self.dL_dK,self.X, self.X2)) def parameters_changed(self): - self.X.gradient[:] = self.kernel.gradients_XX(self.dL_dK, self.X, self.X2,cov=True).sum(0).sum(1) + self.X.gradient[:] = self.kernel.gradients_XX(self.dL_dK, self.X, self.X2) # class Kern_check_d2Kdiag_dXdX(Kern_check_model): # """This class allows gradient checks for the secondderivative of a kernel diagonal with respect to X. """ @@ -263,7 +263,7 @@ def check_kernel_gradient_functions(kern, X=None, X2=None, output_ind=None, verb except NotImplementedError: result=True if verbose: - print(("gradients_XX not implemented for " + kern.name)) + print(("gradients_X not implemented for " + kern.name)) if result and verbose: print("Check passed.") if not result: @@ -283,7 +283,7 @@ def check_kernel_gradient_functions(kern, X=None, X2=None, output_ind=None, verb except NotImplementedError: result=True if verbose: - print(("gradients_XX not implemented for " + kern.name)) + print(("gradients_X not implemented for " + kern.name)) if result and verbose: print("Check passed.") if not result: @@ -293,9 +293,6 @@ def check_kernel_gradient_functions(kern, X=None, X2=None, output_ind=None, verb pass_checks = False return False -# if verbose: -# print("Checking gradients of dKdiag(X, X) wrt X.") - return pass_checks From 0e109cd3dac077156d60fdcfd72420049e27bb06 Mon Sep 17 00:00:00 2001 From: alessandratosi Date: Thu, 21 Apr 2016 15:45:37 +0100 Subject: [PATCH 04/20] syntax fix --- GPy/kern/src/add.py | 5 ++--- GPy/kern/src/kernel_slice_operations.py | 6 +++--- GPy/kern/src/linear.py | 2 +- GPy/kern/src/static.py | 2 +- GPy/kern/src/stationary.py | 8 ++++---- 5 files changed, 11 insertions(+), 12 deletions(-) diff --git a/GPy/kern/src/add.py b/GPy/kern/src/add.py index 9b83d633..7c03d064 100644 --- a/GPy/kern/src/add.py +++ b/GPy/kern/src/add.py @@ -86,7 +86,7 @@ class Add(CombinationKernel): return target def gradients_XX(self, dL_dK, X, X2, cov=True): - if cov==True: # full covarance + if cov: # full covarance if X2 is None: target = np.zeros((X.shape[0], X.shape[0], X.shape[1], X.shape[1])) else: @@ -96,8 +96,7 @@ class Add(CombinationKernel): target = np.zeros((X.shape[0], X.shape[0], X.shape[1])) else: target = np.zeros((X.shape[0], X2.shape[0], X.shape[1])) - - [target.__iadd__(p.gradients_XX(dL_dK, X, X2, cov)) for p in self.parts] + [target.__iadd__(p.gradients_XX(dL_dK, X, X2)) for p in self.parts] return target def gradients_XX_diag(self, dL_dKdiag, X): diff --git a/GPy/kern/src/kernel_slice_operations.py b/GPy/kern/src/kernel_slice_operations.py index ddb16ea1..104bed48 100644 --- a/GPy/kern/src/kernel_slice_operations.py +++ b/GPy/kern/src/kernel_slice_operations.py @@ -119,14 +119,14 @@ def _slice_gradients_XX(f): N, M = X.shape[0], X.shape[0] else: N, M = X.shape[0], X2.shape[0] - if cov==True: # full covariance + if cov: # full covariance with _Slice_wrap(self, X, X2, ret_shape=(N, M, X.shape[1], X.shape[1])) as s: #with _Slice_wrap(self, X, X2, ret_shape=None) as s: - ret = s.handle_return_array(f(self, dL_dK, s.X, s.X2, cov=True)) + ret = s.handle_return_array(f(self, dL_dK, s.X, s.X2, cov)) else: # diagonal covariance with _Slice_wrap(self, X, X2, ret_shape=(N, M, X.shape[1])) as s: #with _Slice_wrap(self, X, X2, ret_shape=None) as s: - ret = s.handle_return_array(f(self, dL_dK, s.X, s.X2, cov=True)) + ret = s.handle_return_array(f(self, dL_dK, s.X, s.X2, cov)) return ret return wrap diff --git a/GPy/kern/src/linear.py b/GPy/kern/src/linear.py index fa412c1d..cd0fb937 100644 --- a/GPy/kern/src/linear.py +++ b/GPy/kern/src/linear.py @@ -101,7 +101,7 @@ class Linear(Kern): #return (((X2[None,:, :] * self.variances)) * dL_dK[:, :, None]).sum(1) return dL_dK.dot(X2)*self.variances #np.einsum('jq,q,ij->iq', X2, self.variances, dL_dK) - def gradients_XX(self, dL_dK, X, X2=None): + def gradients_XX(self, dL_dK, X, X2=None, cov=True): if X2 is None: dL_dK = (dL_dK+dL_dK.T)/2 if X2 is None: return 2*np.ones(X.shape)*self.variances diff --git a/GPy/kern/src/static.py b/GPy/kern/src/static.py index 3ce0dc0a..a56d1903 100644 --- a/GPy/kern/src/static.py +++ b/GPy/kern/src/static.py @@ -24,7 +24,7 @@ class Static(Kern): def gradients_X_diag(self, dL_dKdiag, X): return np.zeros(X.shape) - def gradients_XX(self, dL_dK, X, X2): + def gradients_XX(self, dL_dK, X, X2=None, cov=True): if X2 is None: X2 = X return np.zeros((X.shape[0], X2.shape[0], X.shape[1]), dtype=np.float64) diff --git a/GPy/kern/src/stationary.py b/GPy/kern/src/stationary.py index 4de52d91..ae302266 100644 --- a/GPy/kern/src/stationary.py +++ b/GPy/kern/src/stationary.py @@ -222,14 +222,14 @@ class Stationary(Kern): """ Given the derivative of the objective K(dL_dK), compute the second derivative of K wrt X and X2: - cov = Full: returns the full covariance matrix [QxQ] of the input dimensionfor each pair or vectors - cov = Diag: returns the diagonal of the covariance matrix [QxQ] of the input dimensionfor each pair + cov = True: returns the full covariance matrix [QxQ] of the input dimensionfor each pair or vectors + cov = False: returns the diagonal of the covariance matrix [QxQ] of the input dimensionfor each pair or vectors (computationally more efficient if the full covariance matrix is not needed) ..math: \frac{\partial^2 K}{\partial X2 ^2} = - \frac{\partial^2 K}{\partial X\partial X2} ..returns: - dL2_dXdX2: [NxMxQ] in the cov=Diag case, or [NxMxQxQ] in the cov=full case, + dL2_dXdX2: [NxMxQxQ] in the cov=True case, or [NxMxQ] in the cov=False case, for X [NxQ] and X2[MxQ] (X2 is X if, X2 is None) Thus, we return the second derivative in X2. """ @@ -251,7 +251,7 @@ class Stationary(Kern): # (seems to have a bug: it is subtracted to the first X1 anyway) tmp1[invdist2==0.] -= self.variance - if cov==True: # full covariance + if cov: # full covariance grad = np.empty((X.shape[0], X2.shape[0], X2.shape[1], X.shape[1]), dtype=np.float64) for q in range(self.input_dim): for r in range(self.input_dim): From f7d09f0c759fbf769016ad0c84214baf7e016571 Mon Sep 17 00:00:00 2001 From: alessandratosi Date: Thu, 21 Apr 2016 17:25:10 +0100 Subject: [PATCH 05/20] bug fix --- GPy/kern/src/kernel_slice_operations.py | 6 ++++-- GPy/kern/src/static.py | 5 ++++- GPy/testing/kernel_tests.py | 5 +++-- 3 files changed, 11 insertions(+), 5 deletions(-) diff --git a/GPy/kern/src/kernel_slice_operations.py b/GPy/kern/src/kernel_slice_operations.py index 104bed48..921ac518 100644 --- a/GPy/kern/src/kernel_slice_operations.py +++ b/GPy/kern/src/kernel_slice_operations.py @@ -69,6 +69,8 @@ class _Slice_wrap(object): ret[:, self.k._all_dims_active] = return_val elif len(self.shape) == 3: ret[:, :, self.k._all_dims_active] = return_val + elif len(self.shape) == 4: + ret[:, :, :, self.k._all_dims_active] = return_val return ret return return_val @@ -120,12 +122,12 @@ def _slice_gradients_XX(f): else: N, M = X.shape[0], X2.shape[0] if cov: # full covariance - with _Slice_wrap(self, X, X2, ret_shape=(N, M, X.shape[1], X.shape[1])) as s: #with _Slice_wrap(self, X, X2, ret_shape=None) as s: + with _Slice_wrap(self, X, X2, ret_shape=(N, M, X.shape[1], X.shape[1])) as s: ret = s.handle_return_array(f(self, dL_dK, s.X, s.X2, cov)) else: # diagonal covariance - with _Slice_wrap(self, X, X2, ret_shape=(N, M, X.shape[1])) as s: #with _Slice_wrap(self, X, X2, ret_shape=None) as s: + with _Slice_wrap(self, X, X2, ret_shape=(N, M, X.shape[1])) as s: ret = s.handle_return_array(f(self, dL_dK, s.X, s.X2, cov)) return ret return wrap diff --git a/GPy/kern/src/static.py b/GPy/kern/src/static.py index a56d1903..1745dc23 100644 --- a/GPy/kern/src/static.py +++ b/GPy/kern/src/static.py @@ -27,7 +27,10 @@ class Static(Kern): def gradients_XX(self, dL_dK, X, X2=None, cov=True): if X2 is None: X2 = X - return np.zeros((X.shape[0], X2.shape[0], X.shape[1]), dtype=np.float64) + if cov: + return np.zeros((X.shape[0], X2.shape[0], X.shape[1],X.shape[1]), dtype=np.float64) + else: + return np.zeros((X.shape[0], X2.shape[0], X.shape[1]), dtype=np.float64) def gradients_XX_diag(self, dL_dKdiag, X): return np.zeros(X.shape) diff --git a/GPy/testing/kernel_tests.py b/GPy/testing/kernel_tests.py index a8af50aa..262b7d45 100644 --- a/GPy/testing/kernel_tests.py +++ b/GPy/testing/kernel_tests.py @@ -108,8 +108,8 @@ class Kern_check_d2K_dXdX(Kern_check_model): """This class allows gradient checks for the secondderivative of a kernel with respect to X. """ def __init__(self, kernel=None, dL_dK=None, X=None, X2=None): Kern_check_model.__init__(self,kernel=kernel,dL_dK=dL_dK, X=X, X2=X2) - #self.X = Param('X',X) - #self.link_parameter(self.X) + self.X = Param('X',X) + self.link_parameter(self.X) def log_likelihood(self): return np.sum(self.kernel.gradients_X(self.dL_dK,self.X, self.X2)) @@ -117,6 +117,7 @@ class Kern_check_d2K_dXdX(Kern_check_model): def parameters_changed(self): self.X.gradient[:] = self.kernel.gradients_XX(self.dL_dK, self.X, self.X2) + # class Kern_check_d2Kdiag_dXdX(Kern_check_model): # """This class allows gradient checks for the secondderivative of a kernel diagonal with respect to X. """ From af286ba5280614ecc8371b71ed53cc6447d1183b Mon Sep 17 00:00:00 2001 From: mzwiessele Date: Fri, 22 Apr 2016 15:46:30 +0100 Subject: [PATCH 06/20] [slicing] fixed slicing for second order derivatives --- GPy/core/gp.py | 4 +- GPy/kern/src/add.py | 6 +- GPy/kern/src/kern.py | 22 ++-- GPy/kern/src/kernel_slice_operations.py | 71 +++++++---- GPy/kern/src/linear.py | 14 ++- GPy/kern/src/static.py | 46 ++++--- GPy/kern/src/stationary.py | 12 +- GPy/testing/kernel_tests.py | 157 +++++++++++++++++++++--- 8 files changed, 250 insertions(+), 82 deletions(-) diff --git a/GPy/core/gp.py b/GPy/core/gp.py index 1434573a..1c615cde 100644 --- a/GPy/core/gp.py +++ b/GPy/core/gp.py @@ -378,9 +378,9 @@ class GP(Model): dK_dXnew_full[i] = kern.gradients_X(one, Xnew, self._predictive_variable[[i]]) if full_cov: - dK2_dXdX = kern.gradients_XX(one, Xnew) + dK2_dXdX = kern.gradients_XX(one, Xnew, cov=False) else: - dK2_dXdX = kern.gradients_XX_diag(one, Xnew) + dK2_dXdX = kern.gradients_XX_diag(one, Xnew, cov=False) def compute_cov_inner(wi): if full_cov: diff --git a/GPy/kern/src/add.py b/GPy/kern/src/add.py index 7c03d064..bb04495c 100644 --- a/GPy/kern/src/add.py +++ b/GPy/kern/src/add.py @@ -96,12 +96,12 @@ class Add(CombinationKernel): target = np.zeros((X.shape[0], X.shape[0], X.shape[1])) else: target = np.zeros((X.shape[0], X2.shape[0], X.shape[1])) - [target.__iadd__(p.gradients_XX(dL_dK, X, X2)) for p in self.parts] + [target.__iadd__(p.gradients_XX(dL_dK, X, X2, cov=cov)) for p in self.parts] return target - def gradients_XX_diag(self, dL_dKdiag, X): + def gradients_XX_diag(self, dL_dKdiag, X, cov=True): target = np.zeros(X.shape) - [target.__iadd__(p.gradients_XX_diag(dL_dKdiag, X)) for p in self.parts] + [target.__iadd__(p.gradients_XX_diag(dL_dKdiag, X, cov=cov)) for p in self.parts] return target @Cache_this(limit=3, force_kwargs=['which_parts']) diff --git a/GPy/kern/src/kern.py b/GPy/kern/src/kern.py index 37307e6b..6731a1c3 100644 --- a/GPy/kern/src/kern.py +++ b/GPy/kern/src/kern.py @@ -15,10 +15,10 @@ class Kern(Parameterized): # This adds input slice support. The rather ugly code for slicing can be # found in kernel_slice_operations # __meataclass__ is ignored in Python 3 - needs to be put in the function definiton - #__metaclass__ = KernCallsViaSlicerMeta - #Here, we use the Python module six to support Py3 and Py2 simultaneously + # __metaclass__ = KernCallsViaSlicerMeta + # Here, we use the Python module six to support Py3 and Py2 simultaneously #=========================================================================== - _support_GPU=False + _support_GPU = False def __init__(self, input_dim, active_dims, name, useGPU=False, *a, **kw): """ The base class for a kernel: a positive definite function @@ -62,7 +62,7 @@ class Kern(Parameterized): self.psicomp = PSICOMP_GH() def __setstate__(self, state): - self._all_dims_active = np.arange(0, max(state['active_dims'])+1) + self._all_dims_active = np.arange(0, max(state['active_dims']) + 1) super(Kern, self).__setstate__(state) @property @@ -132,14 +132,14 @@ class Kern(Parameterized): raise NotImplementedError def gradients_X_X2(self, dL_dK, X, X2): return self.gradients_X(dL_dK, X, X2), self.gradients_X(dL_dK.T, X2, X) - def gradients_XX(self, dL_dK, X, X2, cov='True'): + def gradients_XX(self, dL_dK, X, X2, cov=True): """ .. math:: \\frac{\partial^2 L}{\partial X\partial X_2} = \\frac{\partial L}{\partial K}\\frac{\partial^2 K}{\partial X\partial X_2} """ raise(NotImplementedError, "This is the second derivative of K wrt X and X2, and not implemented for this kernel") - def gradients_XX_diag(self, dL_dKdiag, X): + def gradients_XX_diag(self, dL_dKdiag, X, cov=True): """ The diagonal of the second derivative w.r.t. X and X2 """ @@ -292,11 +292,11 @@ class Kern(Parameterized): """ assert isinstance(other, Kern), "only kernels can be multiplied to kernels..." from .prod import Prod - #kernels = [] - #if isinstance(self, Prod): kernels.extend(self.parameters) - #else: kernels.append(self) - #if isinstance(other, Prod): kernels.extend(other.parameters) - #else: kernels.append(other) + # kernels = [] + # if isinstance(self, Prod): kernels.extend(self.parameters) + # else: kernels.append(self) + # if isinstance(other, Prod): kernels.extend(other.parameters) + # else: kernels.append(other) return Prod([self, other], name) def _check_input_dim(self, X): diff --git a/GPy/kern/src/kernel_slice_operations.py b/GPy/kern/src/kernel_slice_operations.py index 921ac518..315f5437 100644 --- a/GPy/kern/src/kernel_slice_operations.py +++ b/GPy/kern/src/kernel_slice_operations.py @@ -25,7 +25,7 @@ class KernCallsViaSlicerMeta(ParametersChangedMeta): put_clean(dct, 'gradients_X', _slice_gradients_X) put_clean(dct, 'gradients_X_X2', _slice_gradients_X) put_clean(dct, 'gradients_XX', _slice_gradients_XX) - put_clean(dct, 'gradients_XX_diag', _slice_gradients_X_diag) + put_clean(dct, 'gradients_XX_diag', _slice_gradients_XX_diag) put_clean(dct, 'gradients_X_diag', _slice_gradients_X_diag) put_clean(dct, 'psi0', _slice_psi) @@ -38,15 +38,16 @@ class KernCallsViaSlicerMeta(ParametersChangedMeta): return super(KernCallsViaSlicerMeta, cls).__new__(cls, name, bases, dct) class _Slice_wrap(object): - def __init__(self, k, X, X2=None, ret_shape=None): + def __init__(self, k, X, X2=None, diag=False, ret_shape=None): self.k = k + self.diag = diag if ret_shape is None: self.shape = X.shape else: self.shape = ret_shape - assert X.ndim == 2, "only matrices are allowed as inputs to kernels for now, given X.shape={!s}".format(X.shape) + assert X.ndim == 2, "need at least column vectors as inputs to kernels for now, given X.shape={!s}".format(X.shape) if X2 is not None: - assert X2.ndim == 2, "only matrices are allowed as inputs to kernels for now, given X2.shape={!s}".format(X2.shape) + assert X2.ndim == 2, "need at least column vectors as inputs to kernels for now, given X2.shape={!s}".format(X2.shape) if (self.k._all_dims_active is not None) and (self.k._sliced_X == 0): self.k._check_active_dims(X) self.X = self.k._slice_X(X) @@ -67,10 +68,13 @@ class _Slice_wrap(object): ret = np.zeros(self.shape) if len(self.shape) == 2: ret[:, self.k._all_dims_active] = return_val - elif len(self.shape) == 3: - ret[:, :, self.k._all_dims_active] = return_val - elif len(self.shape) == 4: - ret[:, :, :, self.k._all_dims_active] = return_val + elif len(self.shape) == 3: # derivative for X2!=None + if self.diag: + ret[:, :, self.k._all_dims_active][:, self.k._all_dims_active] = return_val + else: + ret[:, :, self.k._all_dims_active] = return_val + elif len(self.shape) == 4: # second order derivative + ret[:, :, self.k._all_dims_active][:, :, :, self.k._all_dims_active] = return_val return ret return return_val @@ -114,24 +118,6 @@ def _slice_gradients_X(f): return ret return wrap -def _slice_gradients_XX(f): - @wraps(f) - def wrap(self, dL_dK, X, X2=None, cov=True): - if X2 is None: - N, M = X.shape[0], X.shape[0] - else: - N, M = X.shape[0], X2.shape[0] - if cov: # full covariance - #with _Slice_wrap(self, X, X2, ret_shape=None) as s: - with _Slice_wrap(self, X, X2, ret_shape=(N, M, X.shape[1], X.shape[1])) as s: - ret = s.handle_return_array(f(self, dL_dK, s.X, s.X2, cov)) - else: # diagonal covariance - #with _Slice_wrap(self, X, X2, ret_shape=None) as s: - with _Slice_wrap(self, X, X2, ret_shape=(N, M, X.shape[1])) as s: - ret = s.handle_return_array(f(self, dL_dK, s.X, s.X2, cov)) - return ret - return wrap - def _slice_gradients_X_diag(f): @wraps(f) def wrap(self, dL_dKdiag, X): @@ -140,6 +126,39 @@ def _slice_gradients_X_diag(f): return ret return wrap +def _slice_gradients_XX(f): + @wraps(f) + def wrap(self, dL_dK, X, X2=None, cov=True): + if X2 is None: + N, M = X.shape[0], X.shape[0] + Q1 = Q2 = X.shape[1] + else: + N, M = X.shape[0], X2.shape[0] + Q1, Q2 = X.shape[1], X2.shape[1] + if cov: # full covariance + #with _Slice_wrap(self, X, X2, ret_shape=None) as s: + with _Slice_wrap(self, X, X2, ret_shape=(N, M, Q1, Q2)) as s: + ret = s.handle_return_array(f(self, dL_dK, s.X, s.X2, cov=cov)) + else: # diagonal covariance + #with _Slice_wrap(self, X, X2, ret_shape=None) as s: + with _Slice_wrap(self, X, X2, ret_shape=(N, M, Q1)) as s: + ret = s.handle_return_array(f(self, dL_dK, s.X, s.X2, cov=cov)) + return ret + return wrap + +def _slice_gradients_XX_diag(f): + @wraps(f) + def wrap(self, dL_dKdiag, X, cov=True): + N, Q = X.shape + if cov: # full covariance + with _Slice_wrap(self, X, None, diag=True, ret_shape=(N, Q, Q)) as s: + ret = s.handle_return_array(f(self, dL_dKdiag, s.X, cov=cov)) + else: # diagonal covariance + with _Slice_wrap(self, X, None, ret_shape=(N, Q)) as s: + ret = s.handle_return_array(f(self, dL_dKdiag, s.X, cov=cov)) + return ret + return wrap + def _slice_psi(f): @wraps(f) def wrap(self, Z, variational_posterior): diff --git a/GPy/kern/src/linear.py b/GPy/kern/src/linear.py index cd0fb937..9d9d5933 100644 --- a/GPy/kern/src/linear.py +++ b/GPy/kern/src/linear.py @@ -102,17 +102,21 @@ class Linear(Kern): return dL_dK.dot(X2)*self.variances #np.einsum('jq,q,ij->iq', X2, self.variances, dL_dK) def gradients_XX(self, dL_dK, X, X2=None, cov=True): - if X2 is None: dL_dK = (dL_dK+dL_dK.T)/2 + #if X2 is None: dL_dK = (dL_dK+dL_dK.T)/2 if X2 is None: - return 2*np.ones(X.shape)*self.variances + return 2*self.variances else: - return np.ones(X.shape)*self.variances + return self.variances + def gradients_X_diag(self, dL_dKdiag, X): return 2.*self.variances*dL_dKdiag[:,None]*X - def gradients_XX_diag(self, dL_dKdiag, X): - return 2*np.ones(X.shape)*self.variances + def gradients_XX_diag(self, dL_dKdiag, X, cov=True): + dims = X.shape + if cov: + dims += (X.shape[1],) + return 2*np.ones(dims)*self.variances def input_sensitivity(self, summarize=True): return np.ones(self.input_dim) * self.variances diff --git a/GPy/kern/src/static.py b/GPy/kern/src/static.py index 1745dc23..995f3b5e 100644 --- a/GPy/kern/src/static.py +++ b/GPy/kern/src/static.py @@ -6,6 +6,7 @@ from .kern import Kern import numpy as np from ...core.parameterization import Param from paramz.transformations import Logexp +from paramz.caching import Cache_this class Static(Kern): def __init__(self, input_dim, variance, active_dims, name): @@ -28,11 +29,14 @@ class Static(Kern): if X2 is None: X2 = X if cov: - return np.zeros((X.shape[0], X2.shape[0], X.shape[1],X.shape[1]), dtype=np.float64) + return np.zeros((X.shape[0], X2.shape[0], X.shape[1], X.shape[1]), dtype=np.float64) else: return np.zeros((X.shape[0], X2.shape[0], X.shape[1]), dtype=np.float64) - def gradients_XX_diag(self, dL_dKdiag, X): - return np.zeros(X.shape) + def gradients_XX_diag(self, dL_dKdiag, X, cov=False): + if cov: + return np.zeros((X.shape[0], X.shape[1], X.shape[1]), dtype=np.float64) + else: + return np.zeros(X.shape, dtype=np.float64) def gradients_Z_expectations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior): return np.zeros(Z.shape) @@ -175,17 +179,23 @@ class Fixed(Static): super(Fixed, self).__init__(input_dim, variance, active_dims, name) self.fixed_K = covariance_matrix def K(self, X, X2): - return self.variance * self.fixed_K + if X2 is None: + return self.variance * self.fixed_K + else: + return np.zeros((X.shape[0], X2.shape[0])) def Kdiag(self, X): return self.variance * self.fixed_K.diagonal() def update_gradients_full(self, dL_dK, X, X2=None): - self.variance.gradient = np.einsum('ij,ij', dL_dK, self.fixed_K) + if X2 is None: + self.variance.gradient = np.einsum('ij,ij', dL_dK, self.fixed_K) + else: + self.variance.gradient = 0 def update_gradients_diag(self, dL_dKdiag, X): self.variance.gradient = np.einsum('i,i', dL_dKdiag, np.diagonal(self.fixed_K)) - + def psi2(self, Z, variational_posterior): return np.zeros((Z.shape[0], Z.shape[0]), dtype=np.float64) @@ -227,21 +237,27 @@ class Precomputed(Fixed): :param variance: the variance of the kernel :type variance: float """ + assert input_dim==1, "Precomputed only implemented in one dimension. Use multiple Precomputed kernels to have more dimensions by making use of active_dims" super(Precomputed, self).__init__(input_dim, covariance_matrix, variance, active_dims, name) - def K(self, X, X2=None): + + @Cache_this(limit=2) + def _index(self, X, X2): if X2 is None: - return self.variance * self.fixed_K[X[:,0].astype('int')][:,X[:,0].astype('int')] + i1 = i2 = X.astype('int').flat else: - return self.variance * self.fixed_K[X[:,0].astype('int')][:,X2[:,0].astype('int')] + i1, i2 = X.astype('int').flat, X2.astype('int').flat + return self.fixed_K[i1,:][:,i2] + + def K(self, X, X2=None): + return self.variance * self._index(X, X2) def Kdiag(self, X): - return self.variance * self.fixed_K[X[:,0].astype('int')][:,X[:,0].astype('int')].diagonal() + return self.variance * self._index(X,None).diagonal() def update_gradients_full(self, dL_dK, X, X2=None): - if X2 is None: - self.variance.gradient = np.einsum('ij,ij', dL_dK, self.fixed_K[X[:,0].astype('int')][:,X[:,0].astype('int')]) - else: - self.variance.gradient = np.einsum('ij,ij', dL_dK, self.fixed_K[X[:,0].astype('int')][:,X2[:,0].astype('int')]) + self.variance.gradient = np.einsum('ij,ij', dL_dK, self._index(X, X2)) def update_gradients_diag(self, dL_dKdiag, X): - self.variance.gradient = np.einsum('i,ii', dL_dKdiag, self.fixed_K[X[:,0].astype('int')][:,X[:,0].astype('int')]) \ No newline at end of file + self.variance.gradient = np.einsum('i,ii', dL_dKdiag, self._index(X, None)) + + \ No newline at end of file diff --git a/GPy/kern/src/stationary.py b/GPy/kern/src/stationary.py index ae302266..8f6d1804 100644 --- a/GPy/kern/src/stationary.py +++ b/GPy/kern/src/stationary.py @@ -273,17 +273,19 @@ class Stationary(Kern): #np.sum( - (tmp2*(tmpdist**2)), axis=1, out=grad[:,q]) return grad - def gradients_XX_diag(self, dL_dK, X): + def gradients_XX_diag(self, dL_dK, X, cov=True): """ - Given the derivative of the objective K(dL_dK), compute the second derivative of K wrt X and X2: + Given the derivative of the objective dL_dK, compute the second derivative of K wrt X: ..math: - \frac{\partial^2 K}{\partial X\partial X2} + \frac{\partial^2 K}{\partial X\partial X} ..returns: - dL2_dXdX2: NxMxQ, for X [NxQ] and X2[MxQ] + dL2_dXdX: [NxQ], for X [NxQ] if cov is False, [NxQxQ] if cov is True """ - return np.ones(X.shape) * self.variance/self.lengthscale**2 + if cov: + return np.zeros(X.shape+(X.shape[1],)) + return np.zeros(X.shape)#np.ones(X.shape) * self.variance/self.lengthscale**2 def _gradients_X_pure(self, dL_dK, X, X2=None): invdist = self._inv_dist(X, X2) diff --git a/GPy/testing/kernel_tests.py b/GPy/testing/kernel_tests.py index 262b7d45..f2b95be8 100644 --- a/GPy/testing/kernel_tests.py +++ b/GPy/testing/kernel_tests.py @@ -104,7 +104,7 @@ class Kern_check_dKdiag_dX(Kern_check_dK_dX): def parameters_changed(self): self.X.gradient[:] = self.kernel.gradients_X_diag(self.dL_dK.diagonal(), self.X) -class Kern_check_d2K_dXdX(Kern_check_model): +class Kern_check_d2K_dXdX_cov(Kern_check_model): """This class allows gradient checks for the secondderivative of a kernel with respect to X. """ def __init__(self, kernel=None, dL_dK=None, X=None, X2=None): Kern_check_model.__init__(self,kernel=kernel,dL_dK=dL_dK, X=X, X2=X2) @@ -115,8 +115,55 @@ class Kern_check_d2K_dXdX(Kern_check_model): return np.sum(self.kernel.gradients_X(self.dL_dK,self.X, self.X2)) def parameters_changed(self): - self.X.gradient[:] = self.kernel.gradients_XX(self.dL_dK, self.X, self.X2) + #if self.kernel.name == 'rbf': + # import ipdb;ipdb.set_trace() + grads = self.kernel.gradients_XX(self.dL_dK, self.X, self.X2, cov=True) + self.X.gradient[:] = grads.sum(-1).sum(1) +class Kern_check_d2K_dXdX_no_cov(Kern_check_model): + """This class allows gradient checks for the secondderivative of a kernel with respect to X. """ + def __init__(self, kernel=None, dL_dK=None, X=None, X2=None): + Kern_check_model.__init__(self,kernel=kernel,dL_dK=dL_dK, X=X, X2=X2) + self.X = Param('X',X) + self.link_parameter(self.X) + + def log_likelihood(self): + return np.sum(self.kernel.gradients_X(self.dL_dK,self.X, self.X2)) + + def parameters_changed(self): + #if self.kernel.name == 'rbf': + # import ipdb;ipdb.set_trace() + grads = self.kernel.gradients_XX(self.dL_dK, self.X, self.X2, cov=False) + self.X.gradient[:] = grads.sum(1) + + +class Kern_check_d2Kdiag_dXdX_cov(Kern_check_model): + """This class allows gradient checks for the secondderivative of a kernel with respect to X. """ + def __init__(self, kernel=None, dL_dK=None, X=None, X2=None): + Kern_check_model.__init__(self,kernel=kernel,dL_dK=dL_dK, X=X, X2=X2) + self.X = Param('X',X) + self.link_parameter(self.X) + + def log_likelihood(self): + return np.sum(self.kernel.gradients_X_diag(self.dL_dK.diagonal(),self.X)) + + def parameters_changed(self): + grads = self.kernel.gradients_XX_diag(self.dL_dK.diagonal(), self.X, cov=True) + self.X.gradient[:] = grads.sum(-1) + +class Kern_check_d2Kdiag_dXdX_no_cov(Kern_check_model): + """This class allows gradient checks for the secondderivative of a kernel with respect to X. """ + def __init__(self, kernel=None, dL_dK=None, X=None, X2=None): + Kern_check_model.__init__(self,kernel=kernel,dL_dK=dL_dK, X=X, X2=X2) + self.X = Param('X',X) + self.link_parameter(self.X) + + def log_likelihood(self): + return np.sum(self.kernel.gradients_X_diag(self.dL_dK.diagonal(),self.X)) + + def parameters_changed(self): + grads = self.kernel.gradients_XX_diag(self.dL_dK.diagonal(), self.X, cov=False) + self.X.gradient[:] = grads # class Kern_check_d2Kdiag_dXdX(Kern_check_model): # """This class allows gradient checks for the secondderivative of a kernel diagonal with respect to X. """ @@ -260,7 +307,7 @@ def check_kernel_gradient_functions(kern, X=None, X2=None, output_ind=None, verb if verbose: print("Checking gradients of dK(X, X) wrt X.") try: - testmodel = Kern_check_d2K_dXdX(kern, X=X, X2=None) + testmodel = Kern_check_d2K_dXdX_no_cov(kern, X=X, X2=None) if fixed_X_dims is not None: testmodel.X[:,fixed_X_dims].fix() result = testmodel.checkgrad(verbose=verbose) @@ -276,11 +323,11 @@ def check_kernel_gradient_functions(kern, X=None, X2=None, output_ind=None, verb assert(result) pass_checks = False return False - + if verbose: print("Checking gradients of dK(X, X2) wrt X.") try: - testmodel = Kern_check_d2K_dXdX(kern, X=X, X2=X2) + testmodel = Kern_check_d2K_dXdX_no_cov(kern, X=X, X2=X2) if fixed_X_dims is not None: testmodel.X[:,fixed_X_dims].fix() result = testmodel.checkgrad(verbose=verbose) @@ -297,6 +344,87 @@ def check_kernel_gradient_functions(kern, X=None, X2=None, output_ind=None, verb pass_checks = False return False + if verbose: + print("Checking gradients of dK(X, X) wrt X with full cov in dimensions") + try: + testmodel = Kern_check_d2K_dXdX_cov(kern, X=X, X2=None) + if fixed_X_dims is not None: + testmodel.X[:,fixed_X_dims].fix() + result = testmodel.checkgrad(verbose=verbose) + except NotImplementedError: + result=True + if verbose: + print(("gradients_X not implemented for " + kern.name)) + if result and verbose: + print("Check passed.") + if not result: + print(("Gradient of dK(X, X) wrt X with full cov in dimensions failed for " + kern.name + " covariance function. Gradient values as follows:")) + testmodel.checkgrad(verbose=True) + assert(result) + pass_checks = False + return False + + if verbose: + print("Checking gradients of dK(X, X2) wrt X with full cov in dimensions") + try: + testmodel = Kern_check_d2K_dXdX_cov(kern, X=X, X2=X2) + if fixed_X_dims is not None: + testmodel.X[:,fixed_X_dims].fix() + result = testmodel.checkgrad(verbose=verbose) + except NotImplementedError: + result=True + if verbose: + print(("gradients_X not implemented for " + kern.name)) + if result and verbose: + print("Check passed.") + if not result: + print(("Gradient of dK(X, X2) wrt X failed for " + kern.name + " covariance function. Gradient values as follows:")) + testmodel.checkgrad(verbose=True) + assert(result) + pass_checks = False + return False + + + if verbose: + print("Checking gradients of dKdiag(X, X) wrt X.") + try: + testmodel = Kern_check_d2Kdiag_dXdX_no_cov(kern, X=X, X2=None) + if fixed_X_dims is not None: + testmodel.X[:,fixed_X_dims].fix() + result = testmodel.checkgrad(verbose=verbose) + except NotImplementedError: + result=True + if verbose: + print(("gradients_X not implemented for " + kern.name)) + if result and verbose: + print("Check passed.") + if not result: + print(("Gradient of dKdiag(X, X) wrt X failed for " + kern.name + " covariance function. Gradient values as follows:")) + testmodel.checkgrad(verbose=True) + assert(result) + pass_checks = False + return False + + if verbose: + print("Checking gradients of dKdiag(X, X) wrt X with cov in dimensions") + try: + testmodel = Kern_check_d2Kdiag_dXdX_cov(kern, X=X, X2=None) + if fixed_X_dims is not None: + testmodel.X[:,fixed_X_dims].fix() + result = testmodel.checkgrad(verbose=verbose) + except NotImplementedError: + result=True + if verbose: + print(("gradients_X not implemented for " + kern.name)) + if result and verbose: + print("Check passed.") + if not result: + print(("Gradient of dKdiag(X, X) wrt X with cov in dimensions failed for " + kern.name + " covariance function. Gradient values as follows:")) + testmodel.checkgrad(verbose=True) + assert(result) + pass_checks = False + return False + return pass_checks @@ -304,8 +432,8 @@ def check_kernel_gradient_functions(kern, X=None, X2=None, output_ind=None, verb class KernelGradientTestsContinuous(unittest.TestCase): def setUp(self): self.N, self.D = 10, 5 - self.X = np.random.randn(self.N,self.D) - self.X2 = np.random.randn(self.N+10,self.D) + self.X = np.random.randn(self.N,self.D+1) + self.X2 = np.random.randn(self.N+10,self.D+1) continuous_kerns = ['RBF', 'Linear'] self.kernclasses = [getattr(GPy.kern, s) for s in continuous_kerns] @@ -354,7 +482,7 @@ class KernelGradientTestsContinuous(unittest.TestCase): def test_Add_dims(self): k = GPy.kern.Matern32(2, active_dims=[2,self.D]) + GPy.kern.RBF(2, active_dims=[0,4]) + GPy.kern.Linear(self.D) k.randomize() - self.assertRaises(IndexError, k.K, self.X) + self.assertRaises(IndexError, k.K, self.X[:, :self.D]) k = GPy.kern.Matern32(2, active_dims=[2,self.D-1]) + GPy.kern.RBF(2, active_dims=[0,4]) + GPy.kern.Linear(self.D) k.randomize() # assert it runs: @@ -369,7 +497,7 @@ class KernelGradientTestsContinuous(unittest.TestCase): self.assertTrue(check_kernel_gradient_functions(k, X=self.X, X2=self.X2, verbose=verbose)) def test_RBF(self): - k = GPy.kern.RBF(self.D, ARD=True) + k = GPy.kern.RBF(self.D-1, ARD=True) k.randomize() self.assertTrue(check_kernel_gradient_functions(k, X=self.X, X2=self.X2, verbose=verbose)) @@ -384,9 +512,8 @@ class KernelGradientTestsContinuous(unittest.TestCase): self.assertTrue(check_kernel_gradient_functions(k, X=self.X, X2=self.X2, verbose=verbose)) def test_Fixed(self): - Xall = np.concatenate([self.X, self.X]) - cov = np.dot(Xall, Xall.T) - X = np.arange(self.N).reshape(1,self.N) + cov = np.dot(self.X, self.X.T) + X = np.arange(self.N).reshape(self.N, 1) k = GPy.kern.Fixed(1, cov) k.randomize() self.assertTrue(check_kernel_gradient_functions(k, X=X, X2=None, verbose=verbose)) @@ -409,11 +536,11 @@ class KernelGradientTestsContinuous(unittest.TestCase): def test_Precomputed(self): Xall = np.concatenate([self.X, self.X2]) cov = np.dot(Xall, Xall.T) - X = np.arange(self.N).reshape(1,self.N) - X2 = np.arange(self.N,2*self.N+10).reshape(1,self.N+10) + X = np.arange(self.N).reshape(self.N, 1) + X2 = np.arange(self.N,2*self.N+10).reshape(self.N+10, 1) k = GPy.kern.Precomputed(1, cov) k.randomize() - self.assertTrue(check_kernel_gradient_functions(k, X=X, X2=X2, verbose=verbose)) + self.assertTrue(check_kernel_gradient_functions(k, X=X, X2=X2, verbose=verbose, fixed_X_dims=[0])) class KernelTestsMiscellaneous(unittest.TestCase): def setUp(self): From 401bfbf20c9d065b42fb8e53a1ef720647484828 Mon Sep 17 00:00:00 2001 From: mzwiessele Date: Fri, 22 Apr 2016 16:15:32 +0100 Subject: [PATCH 07/20] [slicing] fixed slicing for second order derivatives --- GPy/kern/src/kern.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/GPy/kern/src/kern.py b/GPy/kern/src/kern.py index 6731a1c3..4379fb71 100644 --- a/GPy/kern/src/kern.py +++ b/GPy/kern/src/kern.py @@ -138,12 +138,12 @@ class Kern(Parameterized): \\frac{\partial^2 L}{\partial X\partial X_2} = \\frac{\partial L}{\partial K}\\frac{\partial^2 K}{\partial X\partial X_2} """ - raise(NotImplementedError, "This is the second derivative of K wrt X and X2, and not implemented for this kernel") + raise NotImplementedError("This is the second derivative of K wrt X and X2, and not implemented for this kernel") def gradients_XX_diag(self, dL_dKdiag, X, cov=True): """ The diagonal of the second derivative w.r.t. X and X2 """ - raise(NotImplementedError, "This is the diagonal of the second derivative of K wrt X and X2, and not implemented for this kernel") + raise NotImplementedError("This is the diagonal of the second derivative of K wrt X and X2, and not implemented for this kernel") def gradients_X_diag(self, dL_dKdiag, X): """ The diagonal of the derivative w.r.t. X From 1cf77c1051db7f45ac27e2b675d469e4887e00f5 Mon Sep 17 00:00:00 2001 From: alessandratosi Date: Mon, 25 Apr 2016 14:53:00 +0100 Subject: [PATCH 08/20] fixed bug in kernel_tests for gradients_XX --- GPy/kern/src/add.py | 2 +- GPy/kern/src/stationary.py | 4 +- GPy/testing/kernel_tests.py | 97 +------------------------------------ 3 files changed, 4 insertions(+), 99 deletions(-) diff --git a/GPy/kern/src/add.py b/GPy/kern/src/add.py index bb04495c..5ac773c9 100644 --- a/GPy/kern/src/add.py +++ b/GPy/kern/src/add.py @@ -100,7 +100,7 @@ class Add(CombinationKernel): return target def gradients_XX_diag(self, dL_dKdiag, X, cov=True): - target = np.zeros(X.shape) + target = np.zeros(X.shape+(X.shape[1],)) [target.__iadd__(p.gradients_XX_diag(dL_dKdiag, X, cov=cov)) for p in self.parts] return target diff --git a/GPy/kern/src/stationary.py b/GPy/kern/src/stationary.py index 8f6d1804..34256d95 100644 --- a/GPy/kern/src/stationary.py +++ b/GPy/kern/src/stationary.py @@ -273,9 +273,9 @@ class Stationary(Kern): #np.sum( - (tmp2*(tmpdist**2)), axis=1, out=grad[:,q]) return grad - def gradients_XX_diag(self, dL_dK, X, cov=True): + def gradients_XX_diag(self, d2L_dK, X, cov=True): """ - Given the derivative of the objective dL_dK, compute the second derivative of K wrt X: + Given the derivative of the objective d2L_dK, compute the second derivative of K wrt X: ..math: \frac{\partial^2 K}{\partial X\partial X} diff --git a/GPy/testing/kernel_tests.py b/GPy/testing/kernel_tests.py index f2b95be8..d89f35e8 100644 --- a/GPy/testing/kernel_tests.py +++ b/GPy/testing/kernel_tests.py @@ -120,25 +120,8 @@ class Kern_check_d2K_dXdX_cov(Kern_check_model): grads = self.kernel.gradients_XX(self.dL_dK, self.X, self.X2, cov=True) self.X.gradient[:] = grads.sum(-1).sum(1) -class Kern_check_d2K_dXdX_no_cov(Kern_check_model): - """This class allows gradient checks for the secondderivative of a kernel with respect to X. """ - def __init__(self, kernel=None, dL_dK=None, X=None, X2=None): - Kern_check_model.__init__(self,kernel=kernel,dL_dK=dL_dK, X=X, X2=X2) - self.X = Param('X',X) - self.link_parameter(self.X) - - def log_likelihood(self): - return np.sum(self.kernel.gradients_X(self.dL_dK,self.X, self.X2)) - - def parameters_changed(self): - #if self.kernel.name == 'rbf': - # import ipdb;ipdb.set_trace() - grads = self.kernel.gradients_XX(self.dL_dK, self.X, self.X2, cov=False) - self.X.gradient[:] = grads.sum(1) - - class Kern_check_d2Kdiag_dXdX_cov(Kern_check_model): - """This class allows gradient checks for the secondderivative of a kernel with respect to X. """ + """This class allows gradient checks for the second derivative of a kernel with respect to X. """ def __init__(self, kernel=None, dL_dK=None, X=None, X2=None): Kern_check_model.__init__(self,kernel=kernel,dL_dK=dL_dK, X=X, X2=X2) self.X = Param('X',X) @@ -151,23 +134,6 @@ class Kern_check_d2Kdiag_dXdX_cov(Kern_check_model): grads = self.kernel.gradients_XX_diag(self.dL_dK.diagonal(), self.X, cov=True) self.X.gradient[:] = grads.sum(-1) -class Kern_check_d2Kdiag_dXdX_no_cov(Kern_check_model): - """This class allows gradient checks for the secondderivative of a kernel with respect to X. """ - def __init__(self, kernel=None, dL_dK=None, X=None, X2=None): - Kern_check_model.__init__(self,kernel=kernel,dL_dK=dL_dK, X=X, X2=X2) - self.X = Param('X',X) - self.link_parameter(self.X) - - def log_likelihood(self): - return np.sum(self.kernel.gradients_X_diag(self.dL_dK.diagonal(),self.X)) - - def parameters_changed(self): - grads = self.kernel.gradients_XX_diag(self.dL_dK.diagonal(), self.X, cov=False) - self.X.gradient[:] = grads - -# class Kern_check_d2Kdiag_dXdX(Kern_check_model): -# """This class allows gradient checks for the secondderivative of a kernel diagonal with respect to X. """ - def check_kernel_gradient_functions(kern, X=None, X2=None, output_ind=None, verbose=False, fixed_X_dims=None): """ This function runs on kernels to check the correctness of their @@ -304,46 +270,6 @@ def check_kernel_gradient_functions(kern, X=None, X2=None, output_ind=None, verb assert(result) return False - if verbose: - print("Checking gradients of dK(X, X) wrt X.") - try: - testmodel = Kern_check_d2K_dXdX_no_cov(kern, X=X, X2=None) - if fixed_X_dims is not None: - testmodel.X[:,fixed_X_dims].fix() - result = testmodel.checkgrad(verbose=verbose) - except NotImplementedError: - result=True - if verbose: - print(("gradients_X not implemented for " + kern.name)) - if result and verbose: - print("Check passed.") - if not result: - print(("Gradient of dK(X, X) wrt X failed for " + kern.name + " covariance function. Gradient values as follows:")) - testmodel.checkgrad(verbose=True) - assert(result) - pass_checks = False - return False - - if verbose: - print("Checking gradients of dK(X, X2) wrt X.") - try: - testmodel = Kern_check_d2K_dXdX_no_cov(kern, X=X, X2=X2) - if fixed_X_dims is not None: - testmodel.X[:,fixed_X_dims].fix() - result = testmodel.checkgrad(verbose=verbose) - except NotImplementedError: - result=True - if verbose: - print(("gradients_X not implemented for " + kern.name)) - if result and verbose: - print("Check passed.") - if not result: - print(("Gradient of dK(X, X2) wrt X failed for " + kern.name + " covariance function. Gradient values as follows:")) - testmodel.checkgrad(verbose=True) - assert(result) - pass_checks = False - return False - if verbose: print("Checking gradients of dK(X, X) wrt X with full cov in dimensions") try: @@ -384,27 +310,6 @@ def check_kernel_gradient_functions(kern, X=None, X2=None, output_ind=None, verb pass_checks = False return False - - if verbose: - print("Checking gradients of dKdiag(X, X) wrt X.") - try: - testmodel = Kern_check_d2Kdiag_dXdX_no_cov(kern, X=X, X2=None) - if fixed_X_dims is not None: - testmodel.X[:,fixed_X_dims].fix() - result = testmodel.checkgrad(verbose=verbose) - except NotImplementedError: - result=True - if verbose: - print(("gradients_X not implemented for " + kern.name)) - if result and verbose: - print("Check passed.") - if not result: - print(("Gradient of dKdiag(X, X) wrt X failed for " + kern.name + " covariance function. Gradient values as follows:")) - testmodel.checkgrad(verbose=True) - assert(result) - pass_checks = False - return False - if verbose: print("Checking gradients of dKdiag(X, X) wrt X with cov in dimensions") try: From bf56a2d85e319fe8aa437497721f3365677aa26f Mon Sep 17 00:00:00 2001 From: alessandratosi Date: Tue, 3 May 2016 15:16:01 +0100 Subject: [PATCH 09/20] fixed covariance computation in predict_jacobian --- GPy/core/gp.py | 37 +++++++++++++++++-------------------- 1 file changed, 17 insertions(+), 20 deletions(-) diff --git a/GPy/core/gp.py b/GPy/core/gp.py index 1c615cde..fc37ab47 100644 --- a/GPy/core/gp.py +++ b/GPy/core/gp.py @@ -335,8 +335,7 @@ class GP(Model): dv_dX += kern.gradients_X(alpha, Xnew, self._predictive_variable) return mean_jac, dv_dX - - def predict_jacobian(self, Xnew, kern=None, full_cov=True): + def predict_jacobian(self, Xnew, kern=None, full_cov=False): """ Compute the derivatives of the posterior of the GP. @@ -354,15 +353,11 @@ class GP(Model): :param X: The points at which to get the predictive gradients. :type X: np.ndarray (Xnew x self.input_dim) :param kern: The kernel to compute the jacobian for. - :param boolean full_cov: whether to return the full covariance of the jacobian. + :param boolean full_cov: whether to return the cross-covariance terms between + the N* Jacobian vectors :returns: dmu_dX, dv_dX :rtype: [np.ndarray (N*, Q ,D), np.ndarray (N*,Q,(D)) ] - - Note: We always return sum in input_dim gradients, as the off-diagonals - in the input_dim are not needed for further calculations. - This is a compromise for increase in speed. Mathematically the jacobian would - have another dimension in Q. """ if kern is None: kern = self.kern @@ -378,27 +373,28 @@ class GP(Model): dK_dXnew_full[i] = kern.gradients_X(one, Xnew, self._predictive_variable[[i]]) if full_cov: - dK2_dXdX = kern.gradients_XX(one, Xnew, cov=False) + dK2_dXdX = kern.gradients_XX(one, Xnew) else: - dK2_dXdX = kern.gradients_XX_diag(one, Xnew, cov=False) + dK2_dXdX = np.zeros((Xnew.shape[0], Xnew.shape[1], Xnew.shape[1])) + for i in range(Xnew.shape[0]): + dK2_dXdX[i:i+1,:,:] = kern.gradients_XX(one, Xnew[i:i+1,:]) def compute_cov_inner(wi): if full_cov: - # full covariance gradients: - var_jac = dK2_dXdX - np.einsum('qnm,miq->niq', dK_dXnew_full.T.dot(wi), dK_dXnew_full) + var_jac = dK2_dXdX - np.einsum('qnm,msr->nsqr', dK_dXnew_full.T.dot(wi), dK_dXnew_full) # n,s = Xnew.shape[0], m = pred_var.shape[0] else: - var_jac = dK2_dXdX - np.einsum('qim,miq->iq', dK_dXnew_full.T.dot(wi), dK_dXnew_full) + var_jac = dK2_dXdX - np.einsum('qnm,mnr->nqr', dK_dXnew_full.T.dot(wi), dK_dXnew_full) return var_jac if self.posterior.woodbury_inv.ndim == 3: # Missing data: if full_cov: - var_jac = np.empty((Xnew.shape[0],Xnew.shape[0],Xnew.shape[1],self.output_dim)) + var_jac = np.empty((Xnew.shape[0],Xnew.shape[0],Xnew.shape[1],Xnew.shape[1],self.output_dim)) + for d in range(self.posterior.woodbury_inv.shape[2]): + var_jac[:, :, :, :, d] = compute_cov_inner(self.posterior.woodbury_inv[:, :, d]) + else: + var_jac = np.empty((Xnew.shape[0],Xnew.shape[1],Xnew.shape[1],self.output_dim)) for d in range(self.posterior.woodbury_inv.shape[2]): var_jac[:, :, :, d] = compute_cov_inner(self.posterior.woodbury_inv[:, :, d]) - else: - var_jac = np.empty((Xnew.shape[0],Xnew.shape[1],self.output_dim)) - for d in range(self.posterior.woodbury_inv.shape[2]): - var_jac[:, :, d] = compute_cov_inner(self.posterior.woodbury_inv[:, :, d]) else: var_jac = compute_cov_inner(self.posterior.woodbury_inv) return mean_jac, var_jac @@ -422,10 +418,11 @@ class GP(Model): mu_jac, var_jac = self.predict_jacobian(Xnew, kern, full_cov=False) mumuT = np.einsum('iqd,ipd->iqp', mu_jac, mu_jac) Sigma = np.zeros(mumuT.shape) - if var_jac.ndim == 3: + if var_jac.ndim == 4: # Missing data Sigma[(slice(None), )+np.diag_indices(Xnew.shape[1], 2)] = var_jac.sum(-1) else: - Sigma[(slice(None), )+np.diag_indices(Xnew.shape[1], 2)] = self.output_dim*var_jac + Sigma = self.output_dim*var_jac + G = 0. if mean: G += mumuT From 5d19039d90de84c4392e741c687a1b7772ca4eb4 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Wed, 4 May 2016 09:11:04 +0100 Subject: [PATCH 10/20] [gradients xx] getting there --- GPy/core/gp.py | 2 +- GPy/kern/src/kernel_slice_operations.py | 4 +- GPy/kern/src/stationary.py | 97 +++++++++++++------------ 3 files changed, 52 insertions(+), 51 deletions(-) diff --git a/GPy/core/gp.py b/GPy/core/gp.py index 1c615cde..dbe66698 100644 --- a/GPy/core/gp.py +++ b/GPy/core/gp.py @@ -449,7 +449,7 @@ class GP(Model): :param bool covariance: whether to include the covariance of the wishart embedding. :param array-like dimensions: which dimensions of the input space to use [defaults to self.get_most_significant_input_dimensions()[:2]] """ - G = self.predict_wishard_embedding(Xnew, kern, mean, covariance) + G = self.predict_wishart_embedding(Xnew, kern, mean, covariance) if dimensions is None: dimensions = self.get_most_significant_input_dimensions()[:2] G = G[:, dimensions][:,:,dimensions] diff --git a/GPy/kern/src/kernel_slice_operations.py b/GPy/kern/src/kernel_slice_operations.py index 315f5437..3f9f6ff3 100644 --- a/GPy/kern/src/kernel_slice_operations.py +++ b/GPy/kern/src/kernel_slice_operations.py @@ -70,11 +70,11 @@ class _Slice_wrap(object): ret[:, self.k._all_dims_active] = return_val elif len(self.shape) == 3: # derivative for X2!=None if self.diag: - ret[:, :, self.k._all_dims_active][:, self.k._all_dims_active] = return_val + ret.T[np.ix_(self.k._all_dims_active, self.k._all_dims_active)] = return_val.T else: ret[:, :, self.k._all_dims_active] = return_val elif len(self.shape) == 4: # second order derivative - ret[:, :, self.k._all_dims_active][:, :, :, self.k._all_dims_active] = return_val + ret.T[np.ix_(self.k._all_dims_active, self.k._all_dims_active)] = return_val.T return ret return return_val diff --git a/GPy/kern/src/stationary.py b/GPy/kern/src/stationary.py index 34256d95..2c93bd68 100644 --- a/GPy/kern/src/stationary.py +++ b/GPy/kern/src/stationary.py @@ -223,7 +223,7 @@ class Stationary(Kern): Given the derivative of the objective K(dL_dK), compute the second derivative of K wrt X and X2: cov = True: returns the full covariance matrix [QxQ] of the input dimensionfor each pair or vectors - cov = False: returns the diagonal of the covariance matrix [QxQ] of the input dimensionfor each pair + cov = False: returns the diagonal of the covariance matrix [QxQ] of the input dimensionfor each pair or vectors (computationally more efficient if the full covariance matrix is not needed) ..math: \frac{\partial^2 K}{\partial X2 ^2} = - \frac{\partial^2 K}{\partial X\partial X2} @@ -247,20 +247,21 @@ class Stationary(Kern): X2 = X tmp1 -= np.eye(X.shape[0])*self.variance else: - #tmp1[X==X2.T] -= self.variance # Old version, to be removed + #tmp1[X==X2.T] -= self.variance # Old version, to be removed # (seems to have a bug: it is subtracted to the first X1 anyway) tmp1[invdist2==0.] -= self.variance - + if cov: # full covariance grad = np.empty((X.shape[0], X2.shape[0], X2.shape[1], X.shape[1]), dtype=np.float64) for q in range(self.input_dim): + tmpdist = (X[:,[q]]-X2[:,[q]].T) for r in range(self.input_dim): - tmpdist2 = (X[:,[q]]-X2[:,[q]].T)*(X[:,[r]]-X2[:,[r]].T) # Introduce temporary distance + tmpdist2 = tmpdist*(X[:,[r]]-X2[:,[r]].T) # Introduce temporary distance if r==q: grad[:, :, q, r] = np.multiply(dL_dK,(np.multiply((tmp1*invdist2 - tmp2),tmpdist2)/l2[r] - tmp1)/l2[q]) else: grad[:, :, q, r] = np.multiply(dL_dK,(np.multiply((tmp1*invdist2 - tmp2),tmpdist2)/l2[r])/l2[q]) - else: + else: # Diagonal covariance, old code grad = np.empty((X.shape[0], X2.shape[0], X.shape[1]), dtype=np.float64) #grad = np.empty(X.shape, dtype=np.float64) @@ -336,18 +337,18 @@ class Exponential(Stationary): def dK_dr(self, r): return -self.K_of_r(r) -# def sde(self): -# """ -# Return the state space representation of the covariance. -# """ -# F = np.array([[-1/self.lengthscale]]) -# L = np.array([[1]]) -# Qc = np.array([[2*self.variance/self.lengthscale]]) -# H = np.array([[1]]) -# Pinf = np.array([[self.variance]]) -# # TODO: return the derivatives as well -# -# return (F, L, Qc, H, Pinf) +# def sde(self): +# """ +# Return the state space representation of the covariance. +# """ +# F = np.array([[-1/self.lengthscale]]) +# L = np.array([[1]]) +# Qc = np.array([[2*self.variance/self.lengthscale]]) +# H = np.array([[1]]) +# Pinf = np.array([[self.variance]]) +# # TODO: return the derivatives as well +# +# return (F, L, Qc, H, Pinf) @@ -416,41 +417,41 @@ class Matern32(Stationary): F1lower = np.array([f(lower) for f in F1])[:, None] return(self.lengthscale ** 3 / (12.*np.sqrt(3) * self.variance) * G + 1. / self.variance * np.dot(Flower, Flower.T) + self.lengthscale ** 2 / (3.*self.variance) * np.dot(F1lower, F1lower.T)) - def sde(self): - """ - Return the state space representation of the covariance. - """ + def sde(self): + """ + Return the state space representation of the covariance. + """ variance = float(self.variance.values) lengthscale = float(self.lengthscale.values) - foo = np.sqrt(3.)/lengthscale - F = np.array([[0, 1], [-foo**2, -2*foo]]) - L = np.array([[0], [1]]) - Qc = np.array([[12.*np.sqrt(3) / lengthscale**3 * variance]]) - H = np.array([[1, 0]]) - Pinf = np.array([[variance, 0], - [0, 3.*variance/(lengthscale**2)]]) - # Allocate space for the derivatives + foo = np.sqrt(3.)/lengthscale + F = np.array([[0, 1], [-foo**2, -2*foo]]) + L = np.array([[0], [1]]) + Qc = np.array([[12.*np.sqrt(3) / lengthscale**3 * variance]]) + H = np.array([[1, 0]]) + Pinf = np.array([[variance, 0], + [0, 3.*variance/(lengthscale**2)]]) + # Allocate space for the derivatives dF = np.empty([F.shape[0],F.shape[1],2]) - dQc = np.empty([Qc.shape[0],Qc.shape[1],2]) - dPinf = np.empty([Pinf.shape[0],Pinf.shape[1],2]) - # The partial derivatives - dFvariance = np.zeros([2,2]) - dFlengthscale = np.array([[0,0], - [6./lengthscale**3,2*np.sqrt(3)/lengthscale**2]]) - dQcvariance = np.array([12.*np.sqrt(3)/lengthscale**3]) - dQclengthscale = np.array([-3*12*np.sqrt(3)/lengthscale**4*variance]) - dPinfvariance = np.array([[1,0],[0,3./lengthscale**2]]) - dPinflengthscale = np.array([[0,0], - [0,-6*variance/lengthscale**3]]) - # Combine the derivatives - dF[:,:,0] = dFvariance - dF[:,:,1] = dFlengthscale - dQc[:,:,0] = dQcvariance - dQc[:,:,1] = dQclengthscale - dPinf[:,:,0] = dPinfvariance - dPinf[:,:,1] = dPinflengthscale + dQc = np.empty([Qc.shape[0],Qc.shape[1],2]) + dPinf = np.empty([Pinf.shape[0],Pinf.shape[1],2]) + # The partial derivatives + dFvariance = np.zeros([2,2]) + dFlengthscale = np.array([[0,0], + [6./lengthscale**3,2*np.sqrt(3)/lengthscale**2]]) + dQcvariance = np.array([12.*np.sqrt(3)/lengthscale**3]) + dQclengthscale = np.array([-3*12*np.sqrt(3)/lengthscale**4*variance]) + dPinfvariance = np.array([[1,0],[0,3./lengthscale**2]]) + dPinflengthscale = np.array([[0,0], + [0,-6*variance/lengthscale**3]]) + # Combine the derivatives + dF[:,:,0] = dFvariance + dF[:,:,1] = dFlengthscale + dQc[:,:,0] = dQcvariance + dQc[:,:,1] = dQclengthscale + dPinf[:,:,0] = dPinfvariance + dPinf[:,:,1] = dPinflengthscale - return (F, L, Qc, H, Pinf, dF, dQc, dPinf) + return (F, L, Qc, H, Pinf, dF, dQc, dPinf) class Matern52(Stationary): """ From e0c71184597f138863cdf9be039c991f79eab651 Mon Sep 17 00:00:00 2001 From: alessandratosi Date: Wed, 4 May 2016 13:07:14 +0100 Subject: [PATCH 11/20] fixed gradients_XX_diag --- GPy/kern/src/stationary.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/GPy/kern/src/stationary.py b/GPy/kern/src/stationary.py index 2c93bd68..519ee80c 100644 --- a/GPy/kern/src/stationary.py +++ b/GPy/kern/src/stationary.py @@ -274,7 +274,7 @@ class Stationary(Kern): #np.sum( - (tmp2*(tmpdist**2)), axis=1, out=grad[:,q]) return grad - def gradients_XX_diag(self, d2L_dK, X, cov=True): + def gradients_XX_diag(self, d2L_dK, X, cov=False): """ Given the derivative of the objective d2L_dK, compute the second derivative of K wrt X: @@ -285,8 +285,9 @@ class Stationary(Kern): dL2_dXdX: [NxQ], for X [NxQ] if cov is False, [NxQxQ] if cov is True """ if cov: - return np.zeros(X.shape+(X.shape[1],)) - return np.zeros(X.shape)#np.ones(X.shape) * self.variance/self.lengthscale**2 + tmp = np.ones(X.shape+(X.shape[1],)) + return tmp * d2L_dK * self.variance/self.lengthscale**2# np.zeros(X.shape+(X.shape[1],)) + return np.ones(X.shape) * d2L_dK * self.variance/self.lengthscale**2 # np.zeros(X.shape) def _gradients_X_pure(self, dL_dK, X, X2=None): invdist = self._inv_dist(X, X2) From b16d57f560fcf5adf68bc6d00bc20899be312d21 Mon Sep 17 00:00:00 2001 From: mzwiessele Date: Thu, 5 May 2016 10:13:46 +0100 Subject: [PATCH 12/20] [dxx] faster numpy version of the gradients_XX --- GPy/kern/src/kernel_slice_operations.py | 2 +- GPy/kern/src/stationary.py | 31 +++++++++++++++---------- 2 files changed, 20 insertions(+), 13 deletions(-) diff --git a/GPy/kern/src/kernel_slice_operations.py b/GPy/kern/src/kernel_slice_operations.py index 3f9f6ff3..73160ef7 100644 --- a/GPy/kern/src/kernel_slice_operations.py +++ b/GPy/kern/src/kernel_slice_operations.py @@ -13,7 +13,7 @@ from paramz.parameterized import ParametersChangedMeta def put_clean(dct, name, func): if name in dct: - #dct['_clean_{}'.format(name)] = dct[name] + dct['_clean_{}'.format(name)] = dct[name] dct[name] = func(dct[name]) class KernCallsViaSlicerMeta(ParametersChangedMeta): diff --git a/GPy/kern/src/stationary.py b/GPy/kern/src/stationary.py index 519ee80c..2a3c1e16 100644 --- a/GPy/kern/src/stationary.py +++ b/GPy/kern/src/stationary.py @@ -237,9 +237,9 @@ class Stationary(Kern): # d2K_dXdX2 = dK_dr*d2r_dXdX2 + d2K_drdr * dr_dX * dr_dX2: invdist = self._inv_dist(X, X2) invdist2 = invdist**2 - dL_dr = self.dK_dr_via_X(X, X2) # * dL_dK we perofrm this product later + dL_dr = self.dK_dr_via_X(X, X2) * dL_dK # we perofrm this product later tmp1 = dL_dr * invdist - dL_drdr = self.dK2_drdr_via_X(X, X2) # * dL_dK we perofrm this product later + dL_drdr = self.dK2_drdr_via_X(X, X2) * dL_dK # we perofrm this product later tmp2 = dL_drdr * invdist2 l2 = np.ones(X.shape[1])*self.lengthscale**2 #np.multiply(np.ones(X.shape[1]) ,self.lengthscale**2) @@ -250,17 +250,24 @@ class Stationary(Kern): #tmp1[X==X2.T] -= self.variance # Old version, to be removed # (seems to have a bug: it is subtracted to the first X1 anyway) tmp1[invdist2==0.] -= self.variance + + tmp3 = (tmp1*invdist2 - tmp2) if cov: # full covariance - grad = np.empty((X.shape[0], X2.shape[0], X2.shape[1], X.shape[1]), dtype=np.float64) - for q in range(self.input_dim): - tmpdist = (X[:,[q]]-X2[:,[q]].T) - for r in range(self.input_dim): - tmpdist2 = tmpdist*(X[:,[r]]-X2[:,[r]].T) # Introduce temporary distance - if r==q: - grad[:, :, q, r] = np.multiply(dL_dK,(np.multiply((tmp1*invdist2 - tmp2),tmpdist2)/l2[r] - tmp1)/l2[q]) - else: - grad[:, :, q, r] = np.multiply(dL_dK,(np.multiply((tmp1*invdist2 - tmp2),tmpdist2)/l2[r])/l2[q]) + dist = X[:,None,:] - X2[None,:,:] + t2 = (tmp3[:,:,None,None]*(dist[:,:,:,None]*dist[:,:,None,:]))/l2[None,None,:,None] + t2.T[np.diag_indices(self.input_dim)] -= tmp1.T[None,:,:] + grad = t2/l2[None,None,None,:] + #grad_old = np.empty((X.shape[0], X2.shape[0], X2.shape[1], X.shape[1]), dtype=np.float64) + # + #for q in range(self.input_dim): + # tmpdist = (X[:,[q]]-X2[:,[q]].T) + # for r in range(self.input_dim): + # tmpdist2 = tmpdist*(X[:,[r]]-X2[:,[r]].T) # Introduce temporary distance + # if r==q: + # grad_old[:, :, q, r] = ((tmp3 * tmpdist2)/l2[r] - tmp1)/l2[q] + # else: + # grad_old[:, :, q, r] = (((tmp3 * tmpdist2)/l2[r])/l2[q]) else: # Diagonal covariance, old code grad = np.empty((X.shape[0], X2.shape[0], X.shape[1]), dtype=np.float64) @@ -274,7 +281,7 @@ class Stationary(Kern): #np.sum( - (tmp2*(tmpdist**2)), axis=1, out=grad[:,q]) return grad - def gradients_XX_diag(self, d2L_dK, X, cov=False): + def gradients_XX_diag(self, d2L_dK, X, cov=True): """ Given the derivative of the objective d2L_dK, compute the second derivative of K wrt X: From 17bfccb45736a1877779218b43791de4e56a3a5e Mon Sep 17 00:00:00 2001 From: mzwiessele Date: Fri, 6 May 2016 16:02:53 +0100 Subject: [PATCH 13/20] [gradxx] not working with X,X... --- GPy/core/gp.py | 2 +- GPy/kern/src/stationary.py | 46 +++++++++++++++++++++++-------------- GPy/plotting/__init__.py | 6 ++++- GPy/testing/kernel_tests.py | 10 ++++---- 4 files changed, 41 insertions(+), 23 deletions(-) diff --git a/GPy/core/gp.py b/GPy/core/gp.py index 78eafa3a..d706ed05 100644 --- a/GPy/core/gp.py +++ b/GPy/core/gp.py @@ -419,7 +419,7 @@ class GP(Model): mumuT = np.einsum('iqd,ipd->iqp', mu_jac, mu_jac) Sigma = np.zeros(mumuT.shape) if var_jac.ndim == 4: # Missing data - Sigma[(slice(None), )+np.diag_indices(Xnew.shape[1], 2)] = var_jac.sum(-1) + Sigma = var_jac.sum(-1) else: Sigma = self.output_dim*var_jac diff --git a/GPy/kern/src/stationary.py b/GPy/kern/src/stationary.py index 2a3c1e16..6ba62fdb 100644 --- a/GPy/kern/src/stationary.py +++ b/GPy/kern/src/stationary.py @@ -237,29 +237,37 @@ class Stationary(Kern): # d2K_dXdX2 = dK_dr*d2r_dXdX2 + d2K_drdr * dr_dX * dr_dX2: invdist = self._inv_dist(X, X2) invdist2 = invdist**2 - dL_dr = self.dK_dr_via_X(X, X2) * dL_dK # we perofrm this product later + dL_dr = self.dK_dr_via_X(X, X2) * dL_dK # we perform this product later tmp1 = dL_dr * invdist dL_drdr = self.dK2_drdr_via_X(X, X2) * dL_dK # we perofrm this product later - tmp2 = dL_drdr * invdist2 + tmp2 = dL_drdr l2 = np.ones(X.shape[1])*self.lengthscale**2 #np.multiply(np.ones(X.shape[1]) ,self.lengthscale**2) - - if X2 is None: - X2 = X - tmp1 -= np.eye(X.shape[0])*self.variance - else: - #tmp1[X==X2.T] -= self.variance # Old version, to be removed - # (seems to have a bug: it is subtracted to the first X1 anyway) - tmp1[invdist2==0.] -= self.variance + + tmp1[invdist2==0.] -= self.variance - tmp3 = (tmp1*invdist2 - tmp2) + tmp3 = (tmp1 - tmp2)*invdist2 + + #tmp3 = (tmp1 - tmp2)*invdist2 + + #tmp3 = tmp3 + # This is not quite right yet, I need the maths to fully understand what is going on.... + #else: if cov: # full covariance - dist = X[:,None,:] - X2[None,:,:] - t2 = (tmp3[:,:,None,None]*(dist[:,:,:,None]*dist[:,:,None,:]))/l2[None,None,:,None] + if X2 is None: + #tmp3 = tmp3+tmp3.T + dist = X[:,None,:] - X[None,:,:] + #dist = dist+dist.swapaxes(0,1) + else: + dist = X[:,None,:] - X2[None,:,:] + dist = (dist[:,:,:,None]*dist[:,:,None,:]) + + t2 = (tmp3[:,:,None,None]*dist)/l2[None,None,:,None] t2.T[np.diag_indices(self.input_dim)] -= tmp1.T[None,:,:] grad = t2/l2[None,None,None,:] + #grad_old = np.empty((X.shape[0], X2.shape[0], X2.shape[1], X.shape[1]), dtype=np.float64) - # + #for q in range(self.input_dim): # tmpdist = (X[:,[q]]-X2[:,[q]].T) # for r in range(self.input_dim): @@ -267,14 +275,18 @@ class Stationary(Kern): # if r==q: # grad_old[:, :, q, r] = ((tmp3 * tmpdist2)/l2[r] - tmp1)/l2[q] # else: - # grad_old[:, :, q, r] = (((tmp3 * tmpdist2)/l2[r])/l2[q]) + # grad_old[:, :, q, r] = (((tmp3 * tmpdist2)/l2[r])/l2[q]) + #import ipdb;ipdb.set_trace() + + if X2 is None: + grad += tmp1[:,:,None,None] else: # Diagonal covariance, old code grad = np.empty((X.shape[0], X2.shape[0], X.shape[1]), dtype=np.float64) #grad = np.empty(X.shape, dtype=np.float64) for q in range(self.input_dim): tmpdist2 = (X[:,[q]]-X2[:,[q]].T) ** 2 - grad[:, :, q] = np.multiply(dL_dK,(np.multiply((tmp1*invdist2 - tmp2),tmpdist2)/l2[q] - tmp1)/l2[q]) + grad[:, :, q] = ((np.multiply((tmp1*invdist2 - tmp2),tmpdist2)/l2[q] - tmp1)/l2[q]) #grad[:, :, q] = ((tmp1*invdist2 - tmp2)*tmpdist2/l2[q] - tmp1)/l2[q] #grad[:, :, q] = ((tmp1*(((tmpdist2)*invdist2/l2[q])-1)) - (tmp2*(tmpdist2))/l2[q])/l2[q] #np.sum(((tmp1*(((tmpdist2)*invdist2/l2[q])-1)) - (tmp2*(tmpdist2))/l2[q])/l2[q], axis=1, out=grad[:,q]) @@ -293,7 +305,7 @@ class Stationary(Kern): """ if cov: tmp = np.ones(X.shape+(X.shape[1],)) - return tmp * d2L_dK * self.variance/self.lengthscale**2# np.zeros(X.shape+(X.shape[1],)) + return tmp * (d2L_dK * self.variance/self.lengthscale**2)[:,None,None]# np.zeros(X.shape+(X.shape[1],)) return np.ones(X.shape) * d2L_dK * self.variance/self.lengthscale**2 # np.zeros(X.shape) def _gradients_X_pure(self, dL_dK, X, X2=None): diff --git a/GPy/plotting/__init__.py b/GPy/plotting/__init__.py index 067f5580..359a841a 100644 --- a/GPy/plotting/__init__.py +++ b/GPy/plotting/__init__.py @@ -50,6 +50,8 @@ def inject_plotting(): GP.plot_samples = gpy_plot.gp_plots.plot_samples GP.plot = gpy_plot.gp_plots.plot GP.plot_f = gpy_plot.gp_plots.plot_f + GP.plot_latent = gpy_plot.gp_plots.plot_f + GP.plot_noiseless = gpy_plot.gp_plots.plot_f GP.plot_magnification = gpy_plot.latent_plots.plot_magnification from ..models import StateSpace @@ -62,7 +64,9 @@ def inject_plotting(): StateSpace.plot_samples = gpy_plot.gp_plots.plot_samples StateSpace.plot = gpy_plot.gp_plots.plot StateSpace.plot_f = gpy_plot.gp_plots.plot_f - + StateSpace.plot_latent = gpy_plot.gp_plots.plot_f + StateSpace.plot_noiseless = gpy_plot.gp_plots.plot_f + from ..core import SparseGP SparseGP.plot_inducing = gpy_plot.data_plots.plot_inducing diff --git a/GPy/testing/kernel_tests.py b/GPy/testing/kernel_tests.py index d89f35e8..b3019de0 100644 --- a/GPy/testing/kernel_tests.py +++ b/GPy/testing/kernel_tests.py @@ -112,13 +112,15 @@ class Kern_check_d2K_dXdX_cov(Kern_check_model): self.link_parameter(self.X) def log_likelihood(self): - return np.sum(self.kernel.gradients_X(self.dL_dK,self.X, self.X2)) + return self.kernel.gradients_X(self.dL_dK, self.X, self.X2).sum() def parameters_changed(self): #if self.kernel.name == 'rbf': - # import ipdb;ipdb.set_trace() - grads = self.kernel.gradients_XX(self.dL_dK, self.X, self.X2, cov=True) - self.X.gradient[:] = grads.sum(-1).sum(1) + # import ipdb;ipdb.set_trace() + if self.X2 is None: X2 = self.X + else: X2 = self.X2 + grads = self.kernel.gradients_XX(self.dL_dK.T, X2, self.X, cov=True) + self.X.gradient[:] = grads.sum(-1).sum(0) class Kern_check_d2Kdiag_dXdX_cov(Kern_check_model): """This class allows gradient checks for the second derivative of a kernel with respect to X. """ From 64f2af719adbd2fe860710c77d256d9154b50db0 Mon Sep 17 00:00:00 2001 From: alessandratosi Date: Fri, 27 May 2016 19:06:28 +0100 Subject: [PATCH 14/20] fixed bug, replaced for loops with einsum --- GPy/kern/src/stationary.py | 54 ++++++++++---------------------------- 1 file changed, 14 insertions(+), 40 deletions(-) diff --git a/GPy/kern/src/stationary.py b/GPy/kern/src/stationary.py index 6ba62fdb..e5093d24 100644 --- a/GPy/kern/src/stationary.py +++ b/GPy/kern/src/stationary.py @@ -237,56 +237,30 @@ class Stationary(Kern): # d2K_dXdX2 = dK_dr*d2r_dXdX2 + d2K_drdr * dr_dX * dr_dX2: invdist = self._inv_dist(X, X2) invdist2 = invdist**2 - dL_dr = self.dK_dr_via_X(X, X2) * dL_dK # we perform this product later + dL_dr = self.dK_dr_via_X(X, X2) #* dL_dK # we perform this product later tmp1 = dL_dr * invdist - dL_drdr = self.dK2_drdr_via_X(X, X2) * dL_dK # we perofrm this product later - tmp2 = dL_drdr + dL_drdr = self.dK2_drdr_via_X(X, X2) #* dL_dK # we perofrm this product later + tmp2 = dL_drdr*invdist2 l2 = np.ones(X.shape[1])*self.lengthscale**2 #np.multiply(np.ones(X.shape[1]) ,self.lengthscale**2) - tmp1[invdist2==0.] -= self.variance - - tmp3 = (tmp1 - tmp2)*invdist2 - - #tmp3 = (tmp1 - tmp2)*invdist2 - - #tmp3 = tmp3 - # This is not quite right yet, I need the maths to fully understand what is going on.... - #else: + if X2 is None: + X2 = X + tmp1 -= np.eye(X.shape[0])*self.variance + else: + tmp1[invdist2==0.] -= self.variance if cov: # full covariance - if X2 is None: - #tmp3 = tmp3+tmp3.T - dist = X[:,None,:] - X[None,:,:] - #dist = dist+dist.swapaxes(0,1) - else: - dist = X[:,None,:] - X2[None,:,:] + grad = np.empty((X.shape[0], X2.shape[0], X2.shape[1], X.shape[1]), dtype=np.float64) + dist = X[:,None,:] - X2[None,:,:] dist = (dist[:,:,:,None]*dist[:,:,None,:]) - - t2 = (tmp3[:,:,None,None]*dist)/l2[None,None,:,None] - t2.T[np.diag_indices(self.input_dim)] -= tmp1.T[None,:,:] - grad = t2/l2[None,None,None,:] - - #grad_old = np.empty((X.shape[0], X2.shape[0], X2.shape[1], X.shape[1]), dtype=np.float64) - - #for q in range(self.input_dim): - # tmpdist = (X[:,[q]]-X2[:,[q]].T) - # for r in range(self.input_dim): - # tmpdist2 = tmpdist*(X[:,[r]]-X2[:,[r]].T) # Introduce temporary distance - # if r==q: - # grad_old[:, :, q, r] = ((tmp3 * tmpdist2)/l2[r] - tmp1)/l2[q] - # else: - # grad_old[:, :, q, r] = (((tmp3 * tmpdist2)/l2[r])/l2[q]) - #import ipdb;ipdb.set_trace() - - if X2 is None: - grad += tmp1[:,:,None,None] - else: - # Diagonal covariance, old code + I = np.ones((X.shape[0], X2.shape[0], X2.shape[1], X.shape[1]))*np.eye((X2.shape[1])) + grad = (np.einsum('kl,klij->klij',dL_dK*(tmp1*invdist2 - tmp2), dist) /l2[None,None,:,None] - np.einsum('kl,klij->klij',dL_dK*tmp1, I))/l2[None,None,None,:] + else: # Diagonal covariance, old code grad = np.empty((X.shape[0], X2.shape[0], X.shape[1]), dtype=np.float64) #grad = np.empty(X.shape, dtype=np.float64) for q in range(self.input_dim): tmpdist2 = (X[:,[q]]-X2[:,[q]].T) ** 2 - grad[:, :, q] = ((np.multiply((tmp1*invdist2 - tmp2),tmpdist2)/l2[q] - tmp1)/l2[q]) + grad[:, :, q] = np.multiply(dL_dK,(np.multiply((tmp1*invdist2 - tmp2),tmpdist2)/l2[q] - tmp1)/l2[q]) #grad[:, :, q] = ((tmp1*invdist2 - tmp2)*tmpdist2/l2[q] - tmp1)/l2[q] #grad[:, :, q] = ((tmp1*(((tmpdist2)*invdist2/l2[q])-1)) - (tmp2*(tmpdist2))/l2[q])/l2[q] #np.sum(((tmp1*(((tmpdist2)*invdist2/l2[q])-1)) - (tmp2*(tmpdist2))/l2[q])/l2[q], axis=1, out=grad[:,q]) From a3f458926b36d06c920a8ea21a49e65113a39baf Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Tue, 7 Jun 2016 09:24:38 +0100 Subject: [PATCH 15/20] [gradsxx] putting tests in, not complete yet! --- GPy/core/gp.py | 11 ++-- GPy/kern/src/kernel_slice_operations.py | 23 +++----- GPy/kern/src/stationary.py | 47 +++++++--------- GPy/testing/kernel_tests.py | 75 +++++++++++++------------ 4 files changed, 72 insertions(+), 84 deletions(-) diff --git a/GPy/core/gp.py b/GPy/core/gp.py index dd9b8b4c..6b1c9cb4 100644 --- a/GPy/core/gp.py +++ b/GPy/core/gp.py @@ -355,7 +355,7 @@ class GP(Model): :param X: The points at which to get the predictive gradients. :type X: np.ndarray (Xnew x self.input_dim) :param kern: The kernel to compute the jacobian for. - :param boolean full_cov: whether to return the cross-covariance terms between + :param boolean full_cov: whether to return the cross-covariance terms between the N* Jacobian vectors :returns: dmu_dX, dv_dX @@ -377,9 +377,10 @@ class GP(Model): if full_cov: dK2_dXdX = kern.gradients_XX(one, Xnew) else: - dK2_dXdX = np.zeros((Xnew.shape[0], Xnew.shape[1], Xnew.shape[1])) - for i in range(Xnew.shape[0]): - dK2_dXdX[i:i+1,:,:] = kern.gradients_XX(one, Xnew[i:i+1,:]) + dK2_dXdX = -kern.gradients_XX(one, Xnew).sum(0) + #dK2_dXdX = np.zeros((Xnew.shape[0], Xnew.shape[1], Xnew.shape[1])) + #for i in range(Xnew.shape[0]): + # dK2_dXdX[i:i+1,:,:] = kern.gradients_XX(one, Xnew[i:i+1,:]) def compute_cov_inner(wi): if full_cov: @@ -424,7 +425,7 @@ class GP(Model): Sigma = var_jac.sum(-1) else: Sigma = self.output_dim*var_jac - + G = 0. if mean: G += mumuT diff --git a/GPy/kern/src/kernel_slice_operations.py b/GPy/kern/src/kernel_slice_operations.py index 73160ef7..c4a64518 100644 --- a/GPy/kern/src/kernel_slice_operations.py +++ b/GPy/kern/src/kernel_slice_operations.py @@ -128,34 +128,25 @@ def _slice_gradients_X_diag(f): def _slice_gradients_XX(f): @wraps(f) - def wrap(self, dL_dK, X, X2=None, cov=True): + def wrap(self, dL_dK, X, X2=None): if X2 is None: N, M = X.shape[0], X.shape[0] Q1 = Q2 = X.shape[1] else: N, M = X.shape[0], X2.shape[0] Q1, Q2 = X.shape[1], X2.shape[1] - if cov: # full covariance - #with _Slice_wrap(self, X, X2, ret_shape=None) as s: - with _Slice_wrap(self, X, X2, ret_shape=(N, M, Q1, Q2)) as s: - ret = s.handle_return_array(f(self, dL_dK, s.X, s.X2, cov=cov)) - else: # diagonal covariance - #with _Slice_wrap(self, X, X2, ret_shape=None) as s: - with _Slice_wrap(self, X, X2, ret_shape=(N, M, Q1)) as s: - ret = s.handle_return_array(f(self, dL_dK, s.X, s.X2, cov=cov)) + #with _Slice_wrap(self, X, X2, ret_shape=None) as s: + with _Slice_wrap(self, X, X2, ret_shape=(N, M, Q1, Q2)) as s: + ret = s.handle_return_array(f(self, dL_dK, s.X, s.X2)) return ret return wrap def _slice_gradients_XX_diag(f): @wraps(f) - def wrap(self, dL_dKdiag, X, cov=True): + def wrap(self, dL_dKdiag, X): N, Q = X.shape - if cov: # full covariance - with _Slice_wrap(self, X, None, diag=True, ret_shape=(N, Q, Q)) as s: - ret = s.handle_return_array(f(self, dL_dKdiag, s.X, cov=cov)) - else: # diagonal covariance - with _Slice_wrap(self, X, None, ret_shape=(N, Q)) as s: - ret = s.handle_return_array(f(self, dL_dKdiag, s.X, cov=cov)) + with _Slice_wrap(self, X, None, diag=True, ret_shape=(N, Q, Q)) as s: + ret = s.handle_return_array(f(self, dL_dKdiag, s.X)) return ret return wrap diff --git a/GPy/kern/src/stationary.py b/GPy/kern/src/stationary.py index e5093d24..22613fa2 100644 --- a/GPy/kern/src/stationary.py +++ b/GPy/kern/src/stationary.py @@ -218,13 +218,13 @@ class Stationary(Kern): else: return self._gradients_X_pure(dL_dK, X, X2) - def gradients_XX(self, dL_dK, X, X2=None, cov=True): + def gradients_XX(self, dL_dK, X, X2=None): """ Given the derivative of the objective K(dL_dK), compute the second derivative of K wrt X and X2: - cov = True: returns the full covariance matrix [QxQ] of the input dimensionfor each pair or vectors - cov = False: returns the diagonal of the covariance matrix [QxQ] of the input dimensionfor each pair - or vectors (computationally more efficient if the full covariance matrix is not needed) + returns the full covariance matrix [QxQ] of the input dimensionfor each pair or vectors, thus + the returned array is of shape [NxNxQxQ]. + ..math: \frac{\partial^2 K}{\partial X2 ^2} = - \frac{\partial^2 K}{\partial X\partial X2} @@ -242,45 +242,36 @@ class Stationary(Kern): dL_drdr = self.dK2_drdr_via_X(X, X2) #* dL_dK # we perofrm this product later tmp2 = dL_drdr*invdist2 l2 = np.ones(X.shape[1])*self.lengthscale**2 #np.multiply(np.ones(X.shape[1]) ,self.lengthscale**2) - + if X2 is None: X2 = X tmp1 -= np.eye(X.shape[0])*self.variance else: tmp1[invdist2==0.] -= self.variance - if cov: # full covariance - grad = np.empty((X.shape[0], X2.shape[0], X2.shape[1], X.shape[1]), dtype=np.float64) - dist = X[:,None,:] - X2[None,:,:] - dist = (dist[:,:,:,None]*dist[:,:,None,:]) - I = np.ones((X.shape[0], X2.shape[0], X2.shape[1], X.shape[1]))*np.eye((X2.shape[1])) - grad = (np.einsum('kl,klij->klij',dL_dK*(tmp1*invdist2 - tmp2), dist) /l2[None,None,:,None] - np.einsum('kl,klij->klij',dL_dK*tmp1, I))/l2[None,None,None,:] - else: # Diagonal covariance, old code - grad = np.empty((X.shape[0], X2.shape[0], X.shape[1]), dtype=np.float64) - #grad = np.empty(X.shape, dtype=np.float64) - for q in range(self.input_dim): - tmpdist2 = (X[:,[q]]-X2[:,[q]].T) ** 2 - grad[:, :, q] = np.multiply(dL_dK,(np.multiply((tmp1*invdist2 - tmp2),tmpdist2)/l2[q] - tmp1)/l2[q]) - #grad[:, :, q] = ((tmp1*invdist2 - tmp2)*tmpdist2/l2[q] - tmp1)/l2[q] - #grad[:, :, q] = ((tmp1*(((tmpdist2)*invdist2/l2[q])-1)) - (tmp2*(tmpdist2))/l2[q])/l2[q] - #np.sum(((tmp1*(((tmpdist2)*invdist2/l2[q])-1)) - (tmp2*(tmpdist2))/l2[q])/l2[q], axis=1, out=grad[:,q]) - #np.sum( - (tmp2*(tmpdist**2)), axis=1, out=grad[:,q]) + #grad = np.empty((X.shape[0], X2.shape[0], X2.shape[1], X.shape[1]), dtype=np.float64) + dist = X[:,None,:] - X2[None,:,:] + dist = (dist[:,:,:,None]*dist[:,:,None,:]) + I = np.ones((X.shape[0], X2.shape[0], X2.shape[1], X.shape[1]))*np.eye((X2.shape[1])) + grad = (np.einsum('kl,klij->klij',dL_dK*(tmp1*invdist2 - tmp2), dist) /l2[None,None,:,None] - np.einsum('kl,klij->klij',dL_dK*tmp1, I))/l2[None,None,None,:] return grad - def gradients_XX_diag(self, d2L_dK, X, cov=True): + def gradients_XX_diag(self, dL_dK_diag, X): """ - Given the derivative of the objective d2L_dK, compute the second derivative of K wrt X: + Given the derivative of the objective dL_dK, compute the second derivative of K wrt X: ..math: \frac{\partial^2 K}{\partial X\partial X} ..returns: - dL2_dXdX: [NxQ], for X [NxQ] if cov is False, [NxQxQ] if cov is True + dL2_dXdX: [NxQxQ] """ - if cov: - tmp = np.ones(X.shape+(X.shape[1],)) - return tmp * (d2L_dK * self.variance/self.lengthscale**2)[:,None,None]# np.zeros(X.shape+(X.shape[1],)) - return np.ones(X.shape) * d2L_dK * self.variance/self.lengthscale**2 # np.zeros(X.shape) + dL_dK_diag = dL_dK_diag.reshape(-1, 1, 1) + assert dL_dK_diag.size == X.shape[0], "dL_dK_diag has to be given as row [N] or column vector [Nx1]" + + l2 = np.ones(X.shape[1])*self.lengthscale**2 + return (dL_dK_diag * self.variance/(l2[:,None]*l2[None,:]))# np.zeros(X.shape+(X.shape[1],)) + #return np.ones(X.shape) * d2L_dK * self.variance/self.lengthscale**2 # np.zeros(X.shape) def _gradients_X_pure(self, dL_dK, X, X2=None): invdist = self._inv_dist(X, X2) diff --git a/GPy/testing/kernel_tests.py b/GPy/testing/kernel_tests.py index b3019de0..9971e6d6 100644 --- a/GPy/testing/kernel_tests.py +++ b/GPy/testing/kernel_tests.py @@ -104,37 +104,42 @@ class Kern_check_dKdiag_dX(Kern_check_dK_dX): def parameters_changed(self): self.X.gradient[:] = self.kernel.gradients_X_diag(self.dL_dK.diagonal(), self.X) -class Kern_check_d2K_dXdX_cov(Kern_check_model): +class Kern_check_d2K_dXdX(Kern_check_model): """This class allows gradient checks for the secondderivative of a kernel with respect to X. """ def __init__(self, kernel=None, dL_dK=None, X=None, X2=None): Kern_check_model.__init__(self,kernel=kernel,dL_dK=dL_dK, X=X, X2=X2) - self.X = Param('X',X) + self.X = Param('X',X.copy()) self.link_parameter(self.X) + self.Xc = X.copy() def log_likelihood(self): + if self.X2 is None: + return self.kernel.gradients_X(self.dL_dK, self.X, self.Xc).sum() return self.kernel.gradients_X(self.dL_dK, self.X, self.X2).sum() def parameters_changed(self): #if self.kernel.name == 'rbf': # import ipdb;ipdb.set_trace() - if self.X2 is None: X2 = self.X - else: X2 = self.X2 - grads = self.kernel.gradients_XX(self.dL_dK.T, X2, self.X, cov=True) - self.X.gradient[:] = grads.sum(-1).sum(0) + if self.X2 is None: + grads = -self.kernel.gradients_XX(self.dL_dK, self.X).sum(1).sum(1) + else: + grads = -self.kernel.gradients_XX(self.dL_dK.T, self.X2, self.X).sum(0).sum(1) + self.X.gradient[:] = grads -class Kern_check_d2Kdiag_dXdX_cov(Kern_check_model): +class Kern_check_d2Kdiag_dXdX(Kern_check_model): """This class allows gradient checks for the second derivative of a kernel with respect to X. """ - def __init__(self, kernel=None, dL_dK=None, X=None, X2=None): - Kern_check_model.__init__(self,kernel=kernel,dL_dK=dL_dK, X=X, X2=X2) + def __init__(self, kernel=None, dL_dK=None, X=None): + Kern_check_model.__init__(self,kernel=kernel,dL_dK=dL_dK, X=X) self.X = Param('X',X) self.link_parameter(self.X) + self.Xc = X.copy() def log_likelihood(self): - return np.sum(self.kernel.gradients_X_diag(self.dL_dK.diagonal(),self.X)) + return np.sum(self.kernel.gradients_X_diag(self.dL_dK.diagonal(), self.X)) def parameters_changed(self): - grads = self.kernel.gradients_XX_diag(self.dL_dK.diagonal(), self.X, cov=True) - self.X.gradient[:] = grads.sum(-1) + grads = self.kernel.gradients_XX_diag(self.dL_dK.diagonal(), self.X) + self.X.gradient[:] = grads.sum(-1) def check_kernel_gradient_functions(kern, X=None, X2=None, output_ind=None, verbose=False, fixed_X_dims=None): """ @@ -273,29 +278,9 @@ def check_kernel_gradient_functions(kern, X=None, X2=None, output_ind=None, verb return False if verbose: - print("Checking gradients of dK(X, X) wrt X with full cov in dimensions") + print("Checking gradients of dK(X, X2) wrt X2 with full cov in dimensions") try: - testmodel = Kern_check_d2K_dXdX_cov(kern, X=X, X2=None) - if fixed_X_dims is not None: - testmodel.X[:,fixed_X_dims].fix() - result = testmodel.checkgrad(verbose=verbose) - except NotImplementedError: - result=True - if verbose: - print(("gradients_X not implemented for " + kern.name)) - if result and verbose: - print("Check passed.") - if not result: - print(("Gradient of dK(X, X) wrt X with full cov in dimensions failed for " + kern.name + " covariance function. Gradient values as follows:")) - testmodel.checkgrad(verbose=True) - assert(result) - pass_checks = False - return False - - if verbose: - print("Checking gradients of dK(X, X2) wrt X with full cov in dimensions") - try: - testmodel = Kern_check_d2K_dXdX_cov(kern, X=X, X2=X2) + testmodel = Kern_check_d2K_dXdX(kern, X=X, X2=X2) if fixed_X_dims is not None: testmodel.X[:,fixed_X_dims].fix() result = testmodel.checkgrad(verbose=verbose) @@ -312,10 +297,30 @@ def check_kernel_gradient_functions(kern, X=None, X2=None, output_ind=None, verb pass_checks = False return False + if verbose: + print("Checking gradients of dK(X, X) wrt X with full cov in dimensions") + try: + testmodel = Kern_check_d2K_dXdX(kern, X=X, X2=None) + if fixed_X_dims is not None: + testmodel.X[:,fixed_X_dims].fix() + result = testmodel.checkgrad(verbose=verbose) + except NotImplementedError: + result=True + if verbose: + print(("gradients_X not implemented for " + kern.name)) + if result and verbose: + print("Check passed.") + if not result: + print(("Gradient of dK(X, X) wrt X with full cov in dimensions failed for " + kern.name + " covariance function. Gradient values as follows:")) + testmodel.checkgrad(verbose=True) + assert(result) + pass_checks = False + return False + if verbose: print("Checking gradients of dKdiag(X, X) wrt X with cov in dimensions") try: - testmodel = Kern_check_d2Kdiag_dXdX_cov(kern, X=X, X2=None) + testmodel = Kern_check_d2Kdiag_dXdX(kern, X=X) if fixed_X_dims is not None: testmodel.X[:,fixed_X_dims].fix() result = testmodel.checkgrad(verbose=verbose) From 787168a3944f42b985c17cd26c6215f447707c8f Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Wed, 8 Jun 2016 10:22:36 +0100 Subject: [PATCH 16/20] [dxxdiag] some steps towards the diagonal gradients in xx --- GPy/core/gp.py | 2 +- GPy/kern/src/add.py | 67 ++++++++++++++++++------------------- GPy/kern/src/stationary.py | 6 ++-- GPy/testing/kernel_tests.py | 7 ++-- 4 files changed, 42 insertions(+), 40 deletions(-) diff --git a/GPy/core/gp.py b/GPy/core/gp.py index 6b1c9cb4..303921ea 100644 --- a/GPy/core/gp.py +++ b/GPy/core/gp.py @@ -377,7 +377,7 @@ class GP(Model): if full_cov: dK2_dXdX = kern.gradients_XX(one, Xnew) else: - dK2_dXdX = -kern.gradients_XX(one, Xnew).sum(0) + dK2_dXdX = kern.gradients_XX(one, Xnew).sum(0) #dK2_dXdX = np.zeros((Xnew.shape[0], Xnew.shape[1], Xnew.shape[1])) #for i in range(Xnew.shape[0]): # dK2_dXdX[i:i+1,:,:] = kern.gradients_XX(one, Xnew[i:i+1,:]) diff --git a/GPy/kern/src/add.py b/GPy/kern/src/add.py index 5ac773c9..99fe809b 100644 --- a/GPy/kern/src/add.py +++ b/GPy/kern/src/add.py @@ -85,23 +85,22 @@ class Add(CombinationKernel): [target.__iadd__(p.gradients_X_diag(dL_dKdiag, X)) for p in self.parts] return target - def gradients_XX(self, dL_dK, X, X2, cov=True): - if cov: # full covarance - if X2 is None: - target = np.zeros((X.shape[0], X.shape[0], X.shape[1], X.shape[1])) - else: - target = np.zeros((X.shape[0], X2.shape[0], X.shape[1], X.shape[1])) - else: # diagonal covariance - if X2 is None: - target = np.zeros((X.shape[0], X.shape[0], X.shape[1])) - else: - target = np.zeros((X.shape[0], X2.shape[0], X.shape[1])) - [target.__iadd__(p.gradients_XX(dL_dK, X, X2, cov=cov)) for p in self.parts] + def gradients_XX(self, dL_dK, X, X2): + if X2 is None: + target = np.zeros((X.shape[0], X.shape[0], X.shape[1], X.shape[1])) + else: + target = np.zeros((X.shape[0], X2.shape[0], X.shape[1], X.shape[1])) + #else: # diagonal covariance + # if X2 is None: + # target = np.zeros((X.shape[0], X.shape[0], X.shape[1])) + # else: + # target = np.zeros((X.shape[0], X2.shape[0], X.shape[1])) + [target.__iadd__(p.gradients_XX(dL_dK, X, X2)) for p in self.parts] return target - def gradients_XX_diag(self, dL_dKdiag, X, cov=True): + def gradients_XX_diag(self, dL_dKdiag, X): target = np.zeros(X.shape+(X.shape[1],)) - [target.__iadd__(p.gradients_XX_diag(dL_dKdiag, X, cov=cov)) for p in self.parts] + [target.__iadd__(p.gradients_XX_diag(dL_dKdiag, X)) for p in self.parts] return target @Cache_this(limit=3, force_kwargs=['which_parts']) @@ -188,7 +187,7 @@ class Add(CombinationKernel): def update_gradients_expectations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior): tmp = dL_dpsi2.sum(0)+ dL_dpsi2.sum(1) if len(dL_dpsi2.shape)==2 else dL_dpsi2.sum(2)+ dL_dpsi2.sum(1) - + if not self._exact_psicomp: return Kern.update_gradients_expectations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior) from .static import White, Bias for p1 in self.parts: @@ -200,9 +199,9 @@ class Add(CombinationKernel): if isinstance(p2, White): continue elif isinstance(p2, Bias): - eff_dL_dpsi1 += tmp * p2.variance + eff_dL_dpsi1 += tmp * p2.variance else:# np.setdiff1d(p1._all_dims_active, ar2, assume_unique): # TODO: Careful, not correct for overlapping _all_dims_active - eff_dL_dpsi1 += tmp * p2.psi1(Z, variational_posterior) + eff_dL_dpsi1 += tmp * p2.psi1(Z, variational_posterior) p1.update_gradients_expectations(dL_dpsi0, eff_dL_dpsi1, dL_dpsi2, Z, variational_posterior) def gradients_Z_expectations(self, dL_psi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior): @@ -219,7 +218,7 @@ class Add(CombinationKernel): if isinstance(p2, White): continue elif isinstance(p2, Bias): - eff_dL_dpsi1 += tmp * p2.variance + eff_dL_dpsi1 += tmp * p2.variance else: eff_dL_dpsi1 += tmp * p2.psi1(Z, variational_posterior) target += p1.gradients_Z_expectations(dL_psi0, eff_dL_dpsi1, dL_dpsi2, Z, variational_posterior) @@ -227,7 +226,7 @@ class Add(CombinationKernel): def gradients_qX_expectations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior): tmp = dL_dpsi2.sum(0)+ dL_dpsi2.sum(1) if len(dL_dpsi2.shape)==2 else dL_dpsi2.sum(2)+ dL_dpsi2.sum(1) - + if not self._exact_psicomp: return Kern.gradients_qX_expectations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior) from .static import White, Bias target_grads = [np.zeros(v.shape) for v in variational_posterior.parameters] @@ -240,9 +239,9 @@ class Add(CombinationKernel): if isinstance(p2, White): continue elif isinstance(p2, Bias): - eff_dL_dpsi1 += tmp * p2.variance + eff_dL_dpsi1 += tmp * p2.variance else: - eff_dL_dpsi1 += tmp * p2.psi1(Z, variational_posterior) + eff_dL_dpsi1 += tmp * p2.psi1(Z, variational_posterior) grads = p1.gradients_qX_expectations(dL_dpsi0, eff_dL_dpsi1, dL_dpsi2, Z, variational_posterior) [np.add(target_grads[i],grads[i],target_grads[i]) for i in range(len(grads))] return target_grads @@ -255,7 +254,7 @@ class Add(CombinationKernel): # other.unlink_parameter(p) # parts.extend(other.parts) # #self.link_parameters(*other_params) - # + # # else: # #self.link_parameter(other) # parts.append(other) @@ -271,7 +270,7 @@ class Add(CombinationKernel): else: return super(Add, self).input_sensitivity(summarize) - + def sde_update_gradient_full(self, gradients): """ Update gradient in the order in which parameters are represented in the @@ -283,12 +282,12 @@ class Add(CombinationKernel): part_param_num = len(p.param_array) # number of parameters in the part p.sde_update_gradient_full(gradients[part_start_param_index:(part_start_param_index+part_param_num)]) part_start_param_index += part_param_num - + def sde(self): """ Support adding kernels for sde representation """ - + import scipy.linalg as la F = None @@ -312,51 +311,51 @@ class Add(CombinationKernel): L = la.block_diag(L,Lt) if (L is not None) else Lt Qc = la.block_diag(Qc,Qct) if (Qc is not None) else Qct H = np.hstack((H,Ht)) if (H is not None) else Ht - + Pinf = la.block_diag(Pinf,Pinft) if (Pinf is not None) else Pinft P0 = la.block_diag(P0,P0t) if (P0 is not None) else P0t - + if dF is not None: dF = np.pad(dF,((0,dFt.shape[0]),(0,dFt.shape[1]),(0,dFt.shape[2])), 'constant', constant_values=0) dF[-dFt.shape[0]:,-dFt.shape[1]:,-dFt.shape[2]:] = dFt else: dF = dFt - + if dQc is not None: dQc = np.pad(dQc,((0,dQct.shape[0]),(0,dQct.shape[1]),(0,dQct.shape[2])), 'constant', constant_values=0) dQc[-dQct.shape[0]:,-dQct.shape[1]:,-dQct.shape[2]:] = dQct else: dQc = dQct - + if dPinf is not None: dPinf = np.pad(dPinf,((0,dPinft.shape[0]),(0,dPinft.shape[1]),(0,dPinft.shape[2])), 'constant', constant_values=0) dPinf[-dPinft.shape[0]:,-dPinft.shape[1]:,-dPinft.shape[2]:] = dPinft else: dPinf = dPinft - + if dP0 is not None: dP0 = np.pad(dP0,((0,dP0t.shape[0]),(0,dP0t.shape[1]),(0,dP0t.shape[2])), 'constant', constant_values=0) dP0[-dP0t.shape[0]:,-dP0t.shape[1]:,-dP0t.shape[2]:] = dP0t else: dP0 = dP0t - + n += Ft.shape[0] nq += Qct.shape[0] nd += dFt.shape[2] - + assert (F.shape[0] == n and F.shape[1]==n), "SDE add: Check of F Dimensions failed" assert (L.shape[0] == n and L.shape[1]==nq), "SDE add: Check of L Dimensions failed" assert (Qc.shape[0] == nq and Qc.shape[1]==nq), "SDE add: Check of Qc Dimensions failed" assert (H.shape[0] == 1 and H.shape[1]==n), "SDE add: Check of H Dimensions failed" assert (Pinf.shape[0] == n and Pinf.shape[1]==n), "SDE add: Check of Pinf Dimensions failed" - assert (P0.shape[0] == n and P0.shape[1]==n), "SDE add: Check of P0 Dimensions failed" + assert (P0.shape[0] == n and P0.shape[1]==n), "SDE add: Check of P0 Dimensions failed" assert (dF.shape[0] == n and dF.shape[1]==n and dF.shape[2]==nd), "SDE add: Check of dF Dimensions failed" assert (dQc.shape[0] == nq and dQc.shape[1]==nq and dQc.shape[2]==nd), "SDE add: Check of dQc Dimensions failed" assert (dPinf.shape[0] == n and dPinf.shape[1]==n and dPinf.shape[2]==nd), "SDE add: Check of dPinf Dimensions failed" assert (dP0.shape[0] == n and dP0.shape[1]==n and dP0.shape[2]==nd), "SDE add: Check of dP0 Dimensions failed" - + return (F,L,Qc,H,Pinf,P0,dF,dQc,dPinf,dP0) diff --git a/GPy/kern/src/stationary.py b/GPy/kern/src/stationary.py index 22613fa2..141a1347 100644 --- a/GPy/kern/src/stationary.py +++ b/GPy/kern/src/stationary.py @@ -266,11 +266,11 @@ class Stationary(Kern): ..returns: dL2_dXdX: [NxQxQ] """ - dL_dK_diag = dL_dK_diag.reshape(-1, 1, 1) + dL_dK_diag = dL_dK_diag.copy().reshape(-1, 1, 1) assert dL_dK_diag.size == X.shape[0], "dL_dK_diag has to be given as row [N] or column vector [Nx1]" - l2 = np.ones(X.shape[1])*self.lengthscale**2 - return (dL_dK_diag * self.variance/(l2[:,None]*l2[None,:]))# np.zeros(X.shape+(X.shape[1],)) + l4 = np.ones(X.shape[1])*self.lengthscale**2 + return dL_dK_diag * (np.eye(X.shape[1]) * self.variance/(l4))[None, :,:]# np.zeros(X.shape+(X.shape[1],)) #return np.ones(X.shape) * d2L_dK * self.variance/self.lengthscale**2 # np.zeros(X.shape) def _gradients_X_pure(self, dL_dK, X, X2=None): diff --git a/GPy/testing/kernel_tests.py b/GPy/testing/kernel_tests.py index 9971e6d6..99951eb1 100644 --- a/GPy/testing/kernel_tests.py +++ b/GPy/testing/kernel_tests.py @@ -135,10 +135,13 @@ class Kern_check_d2Kdiag_dXdX(Kern_check_model): self.Xc = X.copy() def log_likelihood(self): - return np.sum(self.kernel.gradients_X_diag(self.dL_dK.diagonal(), self.X)) + l = 0. + for i in range(self.X.shape[0]): + l += self.kernel.gradients_X(self.dL_dK[[i],[i]], self.X[[i]], self.Xc[[i]]).sum() + return l def parameters_changed(self): - grads = self.kernel.gradients_XX_diag(self.dL_dK.diagonal(), self.X) + grads = -self.kernel.gradients_XX_diag(self.dL_dK.diagonal(), self.X) self.X.gradient[:] = grads.sum(-1) def check_kernel_gradient_functions(kern, X=None, X2=None, output_ind=None, verbose=False, fixed_X_dims=None): From 0c6e3bc88f325280af5bcfa01bc83564afe1b113 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Wed, 8 Jun 2016 13:45:32 +0100 Subject: [PATCH 17/20] [grads x] diagonal entries fixed and add kernel adjusted --- GPy/kern/src/linear.py | 38 ++++++++++++++++++++++++++++---------- GPy/kern/src/rbf.py | 2 ++ GPy/kern/src/static.py | 18 ++++++------------ GPy/kern/src/stationary.py | 10 ++++++++-- 4 files changed, 44 insertions(+), 24 deletions(-) diff --git a/GPy/kern/src/linear.py b/GPy/kern/src/linear.py index 9d9d5933..e7089fe1 100644 --- a/GPy/kern/src/linear.py +++ b/GPy/kern/src/linear.py @@ -101,22 +101,40 @@ class Linear(Kern): #return (((X2[None,:, :] * self.variances)) * dL_dK[:, :, None]).sum(1) return dL_dK.dot(X2)*self.variances #np.einsum('jq,q,ij->iq', X2, self.variances, dL_dK) - def gradients_XX(self, dL_dK, X, X2=None, cov=True): - #if X2 is None: dL_dK = (dL_dK+dL_dK.T)/2 + def gradients_XX(self, dL_dK, X, X2=None): + """ + Given the derivative of the objective K(dL_dK), compute the second derivative of K wrt X and X2: + + returns the full covariance matrix [QxQ] of the input dimensionfor each pair or vectors, thus + the returned array is of shape [NxNxQxQ]. + + ..math: + \frac{\partial^2 K}{\partial X2 ^2} = - \frac{\partial^2 K}{\partial X\partial X2} + + ..returns: + dL2_dXdX2: [NxMxQxQ] for X [NxQ] and X2[MxQ] (X2 is X if, X2 is None) + Thus, we return the second derivative in X2. + """ if X2 is None: - return 2*self.variances - else: - return self.variances + X2 = X + return np.zeros((X.shape[0], X2.shape[0], X.shape[1], X.shape[1])) + #if X2 is None: dL_dK = (dL_dK+dL_dK.T)/2 + #if X2 is None: + # return np.ones(np.repeat(X.shape, 2)) * (self.variances[None,:] + self.variances[:, None])[None, None, :, :] + #else: + # return np.ones((X.shape[0], X2.shape[0], X.shape[1], X.shape[1])) * (self.variances[None,:] + self.variances[:, None])[None, None, :, :] def gradients_X_diag(self, dL_dKdiag, X): return 2.*self.variances*dL_dKdiag[:,None]*X - def gradients_XX_diag(self, dL_dKdiag, X, cov=True): - dims = X.shape - if cov: - dims += (X.shape[1],) - return 2*np.ones(dims)*self.variances + def gradients_XX_diag(self, dL_dKdiag, X): + return np.zeros((X.shape[0], X.shape[1], X.shape[1])) + + #dims = X.shape + #if cov: + # dims += (X.shape[1],) + #return 2*np.ones(dims)*self.variances def input_sensitivity(self, summarize=True): return np.ones(self.input_dim) * self.variances diff --git a/GPy/kern/src/rbf.py b/GPy/kern/src/rbf.py index ff86561d..7a15abe8 100644 --- a/GPy/kern/src/rbf.py +++ b/GPy/kern/src/rbf.py @@ -39,6 +39,8 @@ class RBF(Stationary): def dK2_drdr(self, r): return (r**2-1)*self.K_of_r(r) + def dK2_drdr_diag(self): + return -self.variance # as the diagonal of r is always filled with zeros def __getstate__(self): dc = super(RBF, self).__getstate__() if self.useGPU: diff --git a/GPy/kern/src/static.py b/GPy/kern/src/static.py index 995f3b5e..5cf4a1c9 100644 --- a/GPy/kern/src/static.py +++ b/GPy/kern/src/static.py @@ -25,18 +25,13 @@ class Static(Kern): def gradients_X_diag(self, dL_dKdiag, X): return np.zeros(X.shape) - def gradients_XX(self, dL_dK, X, X2=None, cov=True): + def gradients_XX(self, dL_dK, X, X2=None): if X2 is None: X2 = X - if cov: - return np.zeros((X.shape[0], X2.shape[0], X.shape[1], X.shape[1]), dtype=np.float64) - else: - return np.zeros((X.shape[0], X2.shape[0], X.shape[1]), dtype=np.float64) + return np.zeros((X.shape[0], X2.shape[0], X.shape[1], X.shape[1]), dtype=np.float64) + def gradients_XX_diag(self, dL_dKdiag, X, cov=False): - if cov: - return np.zeros((X.shape[0], X.shape[1], X.shape[1]), dtype=np.float64) - else: - return np.zeros(X.shape, dtype=np.float64) + return np.zeros((X.shape[0], X.shape[1], X.shape[1]), dtype=np.float64) def gradients_Z_expectations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior): return np.zeros(Z.shape) @@ -195,7 +190,7 @@ class Fixed(Static): def update_gradients_diag(self, dL_dKdiag, X): self.variance.gradient = np.einsum('i,i', dL_dKdiag, np.diagonal(self.fixed_K)) - + def psi2(self, Z, variational_posterior): return np.zeros((Z.shape[0], Z.shape[0]), dtype=np.float64) @@ -259,5 +254,4 @@ class Precomputed(Fixed): def update_gradients_diag(self, dL_dKdiag, X): self.variance.gradient = np.einsum('i,ii', dL_dKdiag, self._index(X, None)) - - \ No newline at end of file + diff --git a/GPy/kern/src/stationary.py b/GPy/kern/src/stationary.py index 141a1347..3bf75a4b 100644 --- a/GPy/kern/src/stationary.py +++ b/GPy/kern/src/stationary.py @@ -85,6 +85,11 @@ class Stationary(Kern): def dK2_drdr(self, r): raise NotImplementedError("implement second derivative of covariance wrt r to use this method") + @Cache_this(limit=3, ignore_args=()) + def dK2_drdr_diag(self): + "Second order derivative of K in r_{i,i}. The diagonal entries are always zero, so we do not give it here." + raise NotImplementedError("implement second derivative of covariance wrt r_diag to use this method") + @Cache_this(limit=3, ignore_args=()) def K(self, X, X2=None): """ @@ -253,7 +258,8 @@ class Stationary(Kern): dist = X[:,None,:] - X2[None,:,:] dist = (dist[:,:,:,None]*dist[:,:,None,:]) I = np.ones((X.shape[0], X2.shape[0], X2.shape[1], X.shape[1]))*np.eye((X2.shape[1])) - grad = (np.einsum('kl,klij->klij',dL_dK*(tmp1*invdist2 - tmp2), dist) /l2[None,None,:,None] - np.einsum('kl,klij->klij',dL_dK*tmp1, I))/l2[None,None,None,:] + grad = (((dL_dK*(tmp1*invdist2 - tmp2))[:,:,None,None] * dist)/l2[None,None,:,None] + - (dL_dK*tmp1)[:,:,None,None] * I)/l2[None,None,None,:] return grad def gradients_XX_diag(self, dL_dK_diag, X): @@ -270,7 +276,7 @@ class Stationary(Kern): assert dL_dK_diag.size == X.shape[0], "dL_dK_diag has to be given as row [N] or column vector [Nx1]" l4 = np.ones(X.shape[1])*self.lengthscale**2 - return dL_dK_diag * (np.eye(X.shape[1]) * self.variance/(l4))[None, :,:]# np.zeros(X.shape+(X.shape[1],)) + return dL_dK_diag * (np.eye(X.shape[1]) * -self.dK2_drdr_diag()/(l4))[None, :,:]# np.zeros(X.shape+(X.shape[1],)) #return np.ones(X.shape) * d2L_dK * self.variance/self.lengthscale**2 # np.zeros(X.shape) def _gradients_X_pure(self, dL_dK, X, X2=None): From b1fd7c9aaf5eb05f3e9faf2b415faaca5a08fe04 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Wed, 8 Jun 2016 14:28:25 +0100 Subject: [PATCH 18/20] [grads x] --- GPy/core/gp.py | 2 +- GPy/kern/src/integral.py | 20 ++++++++++---------- GPy/kern/src/stationary.py | 2 +- 3 files changed, 12 insertions(+), 12 deletions(-) diff --git a/GPy/core/gp.py b/GPy/core/gp.py index 303921ea..9cbe3daf 100644 --- a/GPy/core/gp.py +++ b/GPy/core/gp.py @@ -377,7 +377,7 @@ class GP(Model): if full_cov: dK2_dXdX = kern.gradients_XX(one, Xnew) else: - dK2_dXdX = kern.gradients_XX(one, Xnew).sum(0) + dK2_dXdX = kern.gradients_XX_diag(one, Xnew) #dK2_dXdX = np.zeros((Xnew.shape[0], Xnew.shape[1], Xnew.shape[1])) #for i in range(Xnew.shape[0]): # dK2_dXdX[i:i+1,:,:] = kern.gradients_XX(one, Xnew[i:i+1,:]) diff --git a/GPy/kern/src/integral.py b/GPy/kern/src/integral.py index d2827390..971a48a8 100644 --- a/GPy/kern/src/integral.py +++ b/GPy/kern/src/integral.py @@ -10,10 +10,10 @@ class Integral(Kern): #todo do I need to inherit from Stationary """ Integral kernel between... """ - + def __init__(self, input_dim, variances=None, lengthscale=None, ARD=False, active_dims=None, name='integral'): super(Integral, self).__init__(input_dim, active_dims, name) - + if lengthscale is None: lengthscale = np.ones(1) else: @@ -22,7 +22,7 @@ class Integral(Kern): #todo do I need to inherit from Stationary self.lengthscale = Param('lengthscale', lengthscale, Logexp()) #Logexp - transforms to allow positive only values... self.variances = Param('variances', variances, Logexp()) #and here. self.link_parameters(self.variances, self.lengthscale) #this just takes a list of parameters we need to optimise. - + def h(self, z): return 0.5 * z * np.sqrt(math.pi) * math.erf(z) + np.exp(-(z**2)) @@ -36,13 +36,13 @@ class Integral(Kern): #todo do I need to inherit from Stationary for i,x in enumerate(X): for j,x2 in enumerate(X): dK_dl[i,j] = self.variances[0]*self.dk_dl(x[0],x2[0],self.lengthscale[0]) #TODO Multiple length scales - dK_dv[i,j] = self.k_xx(x[0],x2[0],self.lengthscale[0]) #the gradient wrt the variance is k_xx. + dK_dv[i,j] = self.k_xx(x[0],x2[0],self.lengthscale[0]) #the gradient wrt the variance is k_xx. self.lengthscale.gradient = np.sum(dK_dl * dL_dK) self.variances.gradient = np.sum(dK_dv * dL_dK) #print "V%0.5f" % self.variances.gradient #print "L%0.5f" % self.lengthscale.gradient - else: #we're finding dK_xf/Dtheta - print "NEED TO HANDLE TODO!" + else: #we're finding dK_xf/Dtheta + print("NEED TO HANDLE TODO!") #useful little function to help calculate the covariances. def g(self,z): @@ -52,15 +52,15 @@ class Integral(Kern): #todo do I need to inherit from Stationary def k_xx(self,t,tprime,l): return 0.5 * (l**2) * ( self.g(t/l) - self.g((t - tprime)/l) + self.g(tprime/l) - 1) - def k_ff(self,t,tprime,l): + def k_ff(self,t,tprime,l): return np.exp(-((t-tprime)**2)/(l**2)) #rbf - + #covariance between the gradient and the actual value def k_xf(self,t,tprime,l): return 0.5 * np.sqrt(math.pi) * l * (math.erf((t-tprime)/l) + math.erf(tprime/l)) def K(self, X, X2=None): - if X2 is None: + if X2 is None: K_xx = np.zeros([X.shape[0],X.shape[0]]) for i,x in enumerate(X): for j,x2 in enumerate(X): @@ -73,7 +73,7 @@ class Integral(Kern): #todo do I need to inherit from Stationary K_xf[i,j] = self.k_xf(x[0],x2[0],self.lengthscale[0]) #print self.variances[0] return K_xf * self.variances[0] - + def Kdiag(self, X): """I've used the fact that we call this method for K_ff when finding the covariance as a hack so I know if I should return K_ff or K_xx. In this case we're returning K_ff!! diff --git a/GPy/kern/src/stationary.py b/GPy/kern/src/stationary.py index 3bf75a4b..1ce8084f 100644 --- a/GPy/kern/src/stationary.py +++ b/GPy/kern/src/stationary.py @@ -273,7 +273,7 @@ class Stationary(Kern): dL2_dXdX: [NxQxQ] """ dL_dK_diag = dL_dK_diag.copy().reshape(-1, 1, 1) - assert dL_dK_diag.size == X.shape[0], "dL_dK_diag has to be given as row [N] or column vector [Nx1]" + assert (dL_dK_diag.size == X.shape[0]) or (dL_dK_diag.size == 1), "dL_dK_diag has to be given as row [N] or column vector [Nx1]" l4 = np.ones(X.shape[1])*self.lengthscale**2 return dL_dK_diag * (np.eye(X.shape[1]) * -self.dK2_drdr_diag()/(l4))[None, :,:]# np.zeros(X.shape+(X.shape[1],)) From 53169a87872c949a3c90c74d76e85b5537a80cda Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Thu, 9 Jun 2016 08:31:29 +0100 Subject: [PATCH 19/20] [integral] py3 compatability --- GPy/kern/src/integral_limits.py | 30 +++++++++++++++--------------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/GPy/kern/src/integral_limits.py b/GPy/kern/src/integral_limits.py index a1b60859..7006ee6f 100644 --- a/GPy/kern/src/integral_limits.py +++ b/GPy/kern/src/integral_limits.py @@ -10,10 +10,10 @@ class Integral_Limits(Kern): #todo do I need to inherit from Stationary """ Integral kernel, can include limits on each integral value. """ - + def __init__(self, input_dim, variances=None, lengthscale=None, ARD=False, active_dims=None, name='integral'): super(Integral_Limits, self).__init__(input_dim, active_dims, name) - + if lengthscale is None: lengthscale = np.ones(1) else: @@ -22,27 +22,27 @@ class Integral_Limits(Kern): #todo do I need to inherit from Stationary self.lengthscale = Param('lengthscale', lengthscale, Logexp()) #Logexp - transforms to allow positive only values... self.variances = Param('variances', variances, Logexp()) #and here. self.link_parameters(self.variances, self.lengthscale) #this just takes a list of parameters we need to optimise. - + def h(self, z): return 0.5 * z * np.sqrt(math.pi) * math.erf(z) + np.exp(-(z**2)) def dk_dl(self, t, tprime, s, sprime, l): #derivative of the kernel wrt lengthscale return l * ( self.h((t-sprime)/l) - self.h((t - tprime)/l) + self.h((tprime-s)/l) - self.h((s-sprime)/l)) - def update_gradients_full(self, dL_dK, X, X2=None): + def update_gradients_full(self, dL_dK, X, X2=None): if X2 is None: #we're finding dK_xx/dTheta dK_dl = np.zeros([X.shape[0],X.shape[0]]) dK_dv = np.zeros([X.shape[0],X.shape[0]]) for i,x in enumerate(X): for j,x2 in enumerate(X): dK_dl[i,j] = self.variances[0]*self.dk_dl(x[0],x2[0],x[1],x2[1],self.lengthscale[0]) - dK_dv[i,j] = self.k_xx(x[0],x2[0],x[1],x2[1],self.lengthscale[0]) #the gradient wrt the variance is k_xx. + dK_dv[i,j] = self.k_xx(x[0],x2[0],x[1],x2[1],self.lengthscale[0]) #the gradient wrt the variance is k_xx. self.lengthscale.gradient = np.sum(dK_dl * dL_dK) self.variances.gradient = np.sum(dK_dv * dL_dK) #print "V%0.5f" % self.variances.gradient #print "L%0.5f" % self.lengthscale.gradient - else: #we're finding dK_xf/Dtheta - print "NEED TO HANDLE TODO!" + else: #we're finding dK_xf/Dtheta + print("NEED TO HANDLE TODO!") #useful little function to help calculate the covariances. def g(self,z): @@ -50,28 +50,28 @@ class Integral_Limits(Kern): #todo do I need to inherit from Stationary def k_xx(self,t,tprime,s,sprime,l): """Covariance between observed values. - + s and t are one domain of the integral (i.e. the integral between s and t) sprime and tprime are another domain of the integral (i.e. the integral between sprime and tprime) - + We're interested in how correlated these two integrals are. - + Note: We've not multiplied by the variance, this is done in K.""" return 0.5 * (l**2) * ( self.g((t-sprime)/l) + self.g((tprime-s)/l) - self.g((t - tprime)/l) - self.g((s-sprime)/l)) - def k_ff(self,t,tprime,l): + def k_ff(self,t,tprime,l): """Doesn't need s or sprime as we're looking at the 'derivatives', so no domains over which to integrate are required""" return np.exp(-((t-tprime)**2)/(l**2)) #rbf - + def k_xf(self,t,tprime,s,l): """Covariance between the gradient (latent value) and the actual (observed) value. - + Note that sprime isn't actually used in this expression, presumably because the 'primes' are the gradient (latent) values which don't involve an integration, and thus there is no domain over which they're integrated, just a single value that we want.""" return 0.5 * np.sqrt(math.pi) * l * (math.erf((t-tprime)/l) + math.erf((tprime-s)/l)) def K(self, X, X2=None): - if X2 is None: + if X2 is None: K_xx = np.zeros([X.shape[0],X.shape[0]]) for i,x in enumerate(X): for j,x2 in enumerate(X): @@ -83,7 +83,7 @@ class Integral_Limits(Kern): #todo do I need to inherit from Stationary for j,x2 in enumerate(X2): K_xf[i,j] = self.k_xf(x[0],x2[0],x[1],self.lengthscale[0]) #x2[1] unused, see k_xf docstring for explanation. return K_xf * self.variances[0] - + def Kdiag(self, X): """I've used the fact that we call this method for K_ff when finding the covariance as a hack so I know if I should return K_ff or K_xx. In this case we're returning K_ff!! From 1efda71674bfac6a4a494e1ab851d92efe020fb7 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Thu, 9 Jun 2016 10:27:12 +0100 Subject: [PATCH 20/20] [integral] py3 compat --- .../src/multidimensional_integral_limits.py | 46 +++++++++---------- 1 file changed, 23 insertions(+), 23 deletions(-) diff --git a/GPy/kern/src/multidimensional_integral_limits.py b/GPy/kern/src/multidimensional_integral_limits.py index 91983f53..0f473742 100644 --- a/GPy/kern/src/multidimensional_integral_limits.py +++ b/GPy/kern/src/multidimensional_integral_limits.py @@ -10,10 +10,10 @@ class Multidimensional_Integral_Limits(Kern): #todo do I need to inherit from St """ Integral kernel, can include limits on each integral value. """ - + def __init__(self, input_dim, variances=None, lengthscale=None, ARD=False, active_dims=None, name='integral'): super(Multidimensional_Integral_Limits, self).__init__(input_dim, active_dims, name) - + if lengthscale is None: lengthscale = np.ones(1) else: @@ -22,38 +22,38 @@ class Multidimensional_Integral_Limits(Kern): #todo do I need to inherit from St self.lengthscale = Param('lengthscale', lengthscale, Logexp()) #Logexp - transforms to allow positive only values... self.variances = Param('variances', variances, Logexp()) #and here. self.link_parameters(self.variances, self.lengthscale) #this just takes a list of parameters we need to optimise. - + def h(self, z): return 0.5 * z * np.sqrt(math.pi) * math.erf(z) + np.exp(-(z**2)) def dk_dl(self, t, tprime, s, sprime, l): #derivative of the kernel wrt lengthscale return l * ( self.h((t-sprime)/l) - self.h((t - tprime)/l) + self.h((tprime-s)/l) - self.h((s-sprime)/l)) - def update_gradients_full(self, dL_dK, X, X2=None): + def update_gradients_full(self, dL_dK, X, X2=None): #print self.variances if X2 is None: #we're finding dK_xx/dTheta dK_dl_term = np.zeros([X.shape[0],X.shape[0],self.lengthscale.shape[0]]) - k_term = np.zeros([X.shape[0],X.shape[0],self.lengthscale.shape[0]]) - dK_dl = np.zeros([X.shape[0],X.shape[0],self.lengthscale.shape[0]]) + k_term = np.zeros([X.shape[0],X.shape[0],self.lengthscale.shape[0]]) + dK_dl = np.zeros([X.shape[0],X.shape[0],self.lengthscale.shape[0]]) dK_dv = np.zeros([X.shape[0],X.shape[0]]) for il,l in enumerate(self.lengthscale): idx = il*2 for i,x in enumerate(X): for j,x2 in enumerate(X): dK_dl_term[i,j,il] = self.dk_dl(x[idx],x2[idx],x[idx+1],x2[idx+1],l) - k_term[i,j,il] = self.k_xx(x[idx],x2[idx],x[idx+1],x2[idx+1],l) - for il,l in enumerate(self.lengthscale): + k_term[i,j,il] = self.k_xx(x[idx],x2[idx],x[idx+1],x2[idx+1],l) + for il,l in enumerate(self.lengthscale): dK_dl = self.variances[0] * dK_dl_term[:,:,il] for jl, l in enumerate(self.lengthscale): if jl!=il: dK_dl *= k_term[:,:,jl] #dK_dl = np.dot(dK_dl,k_term[:,:,il]) #print k_term[:,:,il] - self.lengthscale.gradient[il] = np.sum(dK_dl * dL_dK) - dK_dv = self.calc_K_xx_wo_variance(X) #the gradient wrt the variance is k_xx. - self.variances.gradient = np.sum(dK_dv * dL_dK) - else: #we're finding dK_xf/Dtheta - print "NEED TO HANDLE TODO!" + self.lengthscale.gradient[il] = np.sum(dK_dl * dL_dK) + dK_dv = self.calc_K_xx_wo_variance(X) #the gradient wrt the variance is k_xx. + self.variances.gradient = np.sum(dK_dv * dL_dK) + else: #we're finding dK_xf/Dtheta + print("NEED TO HANDLE TODO!") #print self.variances[0],self.lengthscale[0],self.lengthscale[1] #np.sum(dK_dv*dL_dK) @@ -63,38 +63,38 @@ class Multidimensional_Integral_Limits(Kern): #todo do I need to inherit from St def k_xx(self,t,tprime,s,sprime,l): """Covariance between observed values. - + s and t are one domain of the integral (i.e. the integral between s and t) sprime and tprime are another domain of the integral (i.e. the integral between sprime and tprime) - + We're interested in how correlated these two integrals are. - + Note: We've not multiplied by the variance, this is done in K.""" return 0.5 * (l**2) * ( self.g((t-sprime)/l) + self.g((tprime-s)/l) - self.g((t - tprime)/l) - self.g((s-sprime)/l)) - def k_ff(self,t,tprime,l): + def k_ff(self,t,tprime,l): """Doesn't need s or sprime as we're looking at the 'derivatives', so no domains over which to integrate are required""" return np.exp(-((t-tprime)**2)/(l**2)) #rbf - + def k_xf(self,t,tprime,s,l): """Covariance between the gradient (latent value) and the actual (observed) value. - + Note that sprime isn't actually used in this expression, presumably because the 'primes' are the gradient (latent) values which don't involve an integration, and thus there is no domain over which they're integrated, just a single value that we want.""" return 0.5 * np.sqrt(math.pi) * l * (math.erf((t-tprime)/l) + math.erf((tprime-s)/l)) def calc_K_xx_wo_variance(self,X): """Calculates K_xx without the variance term""" - K_xx = np.ones([X.shape[0],X.shape[0]]) #ones now as a product occurs over each dimension + K_xx = np.ones([X.shape[0],X.shape[0]]) #ones now as a product occurs over each dimension for i,x in enumerate(X): for j,x2 in enumerate(X): for il,l in enumerate(self.lengthscale): idx = il*2 #each pair of input dimensions describe the limits on one actual dimension in the data K_xx[i,j] *= self.k_xx(x[idx],x2[idx],x[idx+1],x2[idx+1],l) return K_xx - + def K(self, X, X2=None): - if X2 is None: + if X2 is None: #print "X x X" K_xx = self.calc_K_xx_wo_variance(X) return K_xx * self.variances[0] @@ -107,7 +107,7 @@ class Multidimensional_Integral_Limits(Kern): #todo do I need to inherit from St idx = il*2 K_xf[i,j] *= self.k_xf(x[idx],x2[idx],x[idx+1],l) return K_xf * self.variances[0] - + def Kdiag(self, X): """I've used the fact that we call this method for K_ff when finding the covariance as a hack so I know if I should return K_ff or K_xx. In this case we're returning K_ff!!