From 7a8a622b5db1ec45d0dfe3b77b8e97fa1a1c62a1 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Tue, 12 Aug 2014 11:27:50 +0100 Subject: [PATCH 1/4] strange bug in np.einsum fixed when using the _out_ argument (thanks T. Cohn) --- GPy/kern/_src/stationary.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/GPy/kern/_src/stationary.py b/GPy/kern/_src/stationary.py index 5500bb88..5aa570a4 100644 --- a/GPy/kern/_src/stationary.py +++ b/GPy/kern/_src/stationary.py @@ -173,7 +173,7 @@ class Stationary(Kern): tmp *= 2. X2 = X ret = np.empty(X.shape, dtype=np.float64) - [np.einsum('ij,ij->i', tmp, X[:,q][:,None]-X2[:,q][None,:], out=ret[:,q]) for q in xrange(self.input_dim)] + [np.sum(tmp*(X[:,q:q+1]-X2[:,q:q+1]), axis=1, out=ret[:,q]) for q in xrange(self.input_dim)] ret /= self.lengthscale**2 return ret From 3fc603f7780ccfcd714fc6ff4a8fd0c455dc91a8 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Tue, 12 Aug 2014 12:02:16 +0100 Subject: [PATCH 2/4] more bugfix --- GPy/kern/_src/stationary.py | 13 +++---------- 1 file changed, 3 insertions(+), 10 deletions(-) diff --git a/GPy/kern/_src/stationary.py b/GPy/kern/_src/stationary.py index fd7dfd90..67014ab5 100644 --- a/GPy/kern/_src/stationary.py +++ b/GPy/kern/_src/stationary.py @@ -169,17 +169,10 @@ class Stationary(Kern): #the lower memory way with a loop tmp = invdist*dL_dr ret = np.empty(X.shape, dtype=np.float64) -<<<<<<< HEAD - [np.sum(tmp*(X[:,q:q+1]-X2[:,q:q+1]), axis=1, out=ret[:,q]) for q in xrange(self.input_dim)] -======= if X2 is None: - [np.einsum('ij,ij->i', tmp, X[:,q][:,None]-X[:,q][None,:], out=ret[:,q]) for q in xrange(self.input_dim)] - ret2 = np.empty(X.shape, dtype=np.float64) - [np.einsum('ij,ji->j', tmp, X[:,q][:,None]-X[:,q][None,:], out=ret2[:,q]) for q in xrange(self.input_dim)] - ret += ret2 - else: - [np.einsum('ij,ij->i', tmp, X[:,q][:,None]-X2[:,q][None,:], out=ret[:,q]) for q in xrange(self.input_dim)] ->>>>>>> 1061bf52482aa3bf6769db810c955d5fbf51ceae + tmp = tmp + tmp.T + X2 = X + [np.einsum('ij,ij->i', tmp, X[:,q][:,None]-X2[:,q][None,:], out=ret[:,q]) for q in xrange(self.input_dim)] ret /= self.lengthscale**2 return ret From fe405271e5a498e118148629b399c1a11646f83a Mon Sep 17 00:00:00 2001 From: James Hensman Date: Wed, 13 Aug 2014 10:36:54 +0100 Subject: [PATCH 3/4] gradients of predictions for Trevor --- GPy/core/gp.py | 22 ++++++++++++++++++++++ GPy/kern/_src/stationary.py | 18 +++++++++--------- GPy/models/gplvm.py | 6 +++--- 3 files changed, 34 insertions(+), 12 deletions(-) diff --git a/GPy/core/gp.py b/GPy/core/gp.py index e0f5755c..ba7c38f2 100644 --- a/GPy/core/gp.py +++ b/GPy/core/gp.py @@ -148,6 +148,28 @@ class GP(Model): m, v = self._raw_predict(X, full_cov=False) return self.likelihood.predictive_quantiles(m, v, quantiles, Y_metadata) + def predictive_gradients(self, Xnew): + """ + Compute the derivatives of the latent function with respect to X* + + Given a set of points at which to predict X* (size [N*,Q]), compute the + derivatives of the mean and variance. Resulting arrays are sized: + dmu_dX* -- [N*, Q ,D], where D is the number of output in this GP (usually one). + dv_dX* -- [N*, Q], (since all outputs have the same variance) + + """ + dmu_dX = np.empty((Xnew.shape[0],Xnew.shape[1],self.output_dim)) + for i in range(self.output_dim): + dmu_dX[:,:,i] = self.kern.gradients_X(self.posterior.woodbury_vector[:,i:i+1].T, Xnew, self.X) + + # gradients wrt the diagonal part k_{xx} + dv_dX = self.kern.gradients_X(np.eye(Xnew.shape[0]), Xnew) + #grads wrt 'Schur' part K_{xf}K_{ff}^{-1}K_{fx} + alpha = -2.*np.dot(self.kern.K(Xnew, self.X),self.posterior.woodbury_inv) + dv_dX += self.kern.gradients_X(alpha, Xnew, self.X) + return dmu_dX, dv_dX + + def posterior_samples_f(self,X,size=10, full_cov=True): """ Samples the posterior GP at the points X. diff --git a/GPy/kern/_src/stationary.py b/GPy/kern/_src/stationary.py index 67014ab5..c0553238 100644 --- a/GPy/kern/_src/stationary.py +++ b/GPy/kern/_src/stationary.py @@ -160,20 +160,20 @@ class Stationary(Kern): """ invdist = self._inv_dist(X, X2) dL_dr = self.dK_dr_via_X(X, X2) * dL_dK - #The high-memory numpy way: - #d = X[:, None, :] - X2[None, :, :] - #ret = np.sum((invdist*dL_dr)[:,:,None]*d,1)/self.lengthscale**2 - #if X2 is None: - #ret *= 2. - - #the lower memory way with a loop tmp = invdist*dL_dr - ret = np.empty(X.shape, dtype=np.float64) if X2 is None: tmp = tmp + tmp.T X2 = X - [np.einsum('ij,ij->i', tmp, X[:,q][:,None]-X2[:,q][None,:], out=ret[:,q]) for q in xrange(self.input_dim)] + + #The high-memory numpy way: + #d = X[:, None, :] - X2[None, :, :] + #ret = np.sum(tmp[:,:,None]*d,1)/self.lengthscale**2 + + #the lower memory way with a loop + ret = np.empty(X.shape, dtype=np.float64) + [np.sum(tmp*(X[:,q][:,None]-X2[:,q][None,:]), axis=1, out=ret[:,q]) for q in xrange(self.input_dim)] ret /= self.lengthscale**2 + return ret def gradients_X_diag(self, dL_dKdiag, X): diff --git a/GPy/models/gplvm.py b/GPy/models/gplvm.py index 542dcd31..8f5432ba 100644 --- a/GPy/models/gplvm.py +++ b/GPy/models/gplvm.py @@ -45,10 +45,10 @@ class GPLVM(GP): self.X.gradient = self.kern.gradients_X(self.grad_dict['dL_dK'], self.X, None) def jacobian(self,X): - target = np.zeros((X.shape[0],X.shape[1],self.output_dim)) + J = np.zeros((X.shape[0],X.shape[1],self.output_dim)) for i in range(self.output_dim): - target[:,:,i]=self.kern.gradients_X(np.dot(self.Ki,self.likelihood.Y[:,i])[None, :],X,self.X) - return target + J[:,:,i] = self.kern.gradients_X(self.posterior.woodbury_vector[:,i:i+1], X, self.X) + return J def magnification(self,X): target=np.zeros(X.shape[0]) From 1b9995776b9a74cb43d622e6f0560078a137ec66 Mon Sep 17 00:00:00 2001 From: Ricardo Date: Wed, 13 Aug 2014 13:51:39 +0100 Subject: [PATCH 4/4] minor changes in ICM --- GPy/util/multioutput.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/GPy/util/multioutput.py b/GPy/util/multioutput.py index b769d6a5..cc9af29e 100644 --- a/GPy/util/multioutput.py +++ b/GPy/util/multioutput.py @@ -40,7 +40,7 @@ def build_likelihood(Y_list,noise_index,likelihoods_list=None): return likelihood -def ICM(input_dim, num_outputs, kernel, W_rank=1,W=None,kappa=None,name='X'): +def ICM(input_dim, num_outputs, kernel, W_rank=1,W=None,kappa=None,name='ICM'): """ Builds a kernel for an Intrinsic Coregionalization Model @@ -56,12 +56,10 @@ def ICM(input_dim, num_outputs, kernel, W_rank=1,W=None,kappa=None,name='X'): warnings.warn("kernel's input dimension overwritten to fit input_dim parameter.") K = kernel.prod(GPy.kern.Coregionalize(1, num_outputs, active_dims=[input_dim], rank=W_rank,W=W,kappa=kappa,name='B'),name=name) - K['.*variance'] = 1. - K['.*variance'].fix() return K -def LCM(input_dim, num_outputs, kernels_list, W_rank=1,name='X'): +def LCM(input_dim, num_outputs, kernels_list, W_rank=1,name='ICM'): """ Builds a kernel for an Linear Coregionalization Model @@ -77,6 +75,7 @@ def LCM(input_dim, num_outputs, kernels_list, W_rank=1,name='X'): j = 1 for kernel in kernels_list[1:]: K += ICM(input_dim,num_outputs,kernel,W_rank,name='%s%s' %(name,j)) + j += 1 return K