From 8b73dafbaeeed497504e0075dc231d1ea2fe9122 Mon Sep 17 00:00:00 2001 From: Nicolo Fusi Date: Mon, 3 Dec 2012 10:08:12 +0000 Subject: [PATCH 1/3] added precomputation of linear kernel, changed the logic a bit --- GPy/kern/linear.py | 22 ++++++++++++++++++++-- 1 file changed, 20 insertions(+), 2 deletions(-) diff --git a/GPy/kern/linear.py b/GPy/kern/linear.py index d3e3f42c..7246244e 100644 --- a/GPy/kern/linear.py +++ b/GPy/kern/linear.py @@ -21,6 +21,7 @@ class linear(kernpart): self.Nparam = 1 self.name = 'linear' self.set_param(variance) + self._Xcache, self._X2cache = np.empty(shape=(2,)) def get_param(self): return self.variance @@ -32,7 +33,8 @@ class linear(kernpart): return ['variance'] def K(self,X,X2,target): - target += self.variance * np.dot(X, X2.T) + self._K_computations(X, X2) + target += self.variance * self._dot_product def Kdiag(self,X,target): np.add(target,np.sum(self.variance*np.square(X),-1),target) @@ -42,7 +44,9 @@ class linear(kernpart): Computes the derivatives wrt theta Return shape is NxMx(Ntheta) """ - product = np.dot(X, X2.T) + self._K_computations(X, X2) + product = self._dot_product + # product = np.dot(X, X2.T) target += np.sum(product*partial) def dK_dX(self,partial,X,X2,target): @@ -51,6 +55,20 @@ class linear(kernpart): def dKdiag_dtheta(self,partial,X,target): target += np.sum(partial*np.square(X).sum(1)) + def _K_computations(self,X,X2): + # (Nicolo) changed the logic here. If X2 is None, we want to cache + # (X,X). In practice X2 should always be passed. + if X2 is None: + X2 = X + if not (np.all(X==self._Xcache) and np.all(X2==self._X2cache)): + self._Xcache = X + self._X2cache = X2 + self._dot_product = np.dot(X,X2.T) + else: + # print "Cache hit!" + pass # TODO: insert debug message here (logging framework) + + # def psi0(self,Z,mu,S,target): # expected = np.square(mu) + S # np.add(target,np.sum(self.variance*expected),target) From 3edd867ece4fa503d537e5488581718538429b8f Mon Sep 17 00:00:00 2001 From: Nicolas Date: Mon, 3 Dec 2012 10:26:28 +0000 Subject: [PATCH 2/3] GPy: Some rewriting for the exponential and Matern kernels. They now pass the unit test. --- GPy/kern/Matern32.py | 51 +++++++++++++++++++++----------------- GPy/kern/Matern52.py | 55 +++++++++++++++++++++++------------------ GPy/kern/exponential.py | 6 ++--- 3 files changed, 63 insertions(+), 49 deletions(-) diff --git a/GPy/kern/Matern32.py b/GPy/kern/Matern32.py index c37bd5c0..8223b37a 100644 --- a/GPy/kern/Matern32.py +++ b/GPy/kern/Matern32.py @@ -39,11 +39,13 @@ class Matern32(kernpart): def get_param(self): """return the value of the parameters.""" return np.hstack((self.variance,self.lengthscales)) + def set_param(self,x): """set the value of the parameters.""" assert x.size==(self.D+1) self.variance = x[0] self.lengthscales = x[1:] + def get_param_names(self): """return parameter names.""" if self.D==1: @@ -56,10 +58,37 @@ class Matern32(kernpart): if X2 is None: X2 = X dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscales),-1)) np.add(self.variance*(1+np.sqrt(3.)*dist)*np.exp(-np.sqrt(3.)*dist), target,target) + def Kdiag(self,X,target): """Compute the diagonal of the covariance matrix associated to X.""" np.add(target,self.variance,target) + def dK_dtheta(self,partial,X,X2,target): + """derivative of the covariance matrix with respect to the parameters.""" + if X2 is None: X2 = X + dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscales),-1)) + dvar = (1+np.sqrt(3.)*dist)*np.exp(-np.sqrt(3.)*dist) + invdist = 1./np.where(dist!=0.,dist,np.inf) + dist2M = np.square(X[:,None,:]-X2[None,:,:])/self.lengthscales**3 + dl = (self.variance* 3 * dist * np.exp(-np.sqrt(3.)*dist))[:,:,np.newaxis] * dist2M*invdist[:,:,np.newaxis] + target[0] += np.sum(dvar*partial) + target[1:] += (dl*partial[:,:,None]).sum(0).sum(0) + + def dKdiag_dtheta(self,partial,X,target): + """derivative of the diagonal of the covariance matrix with respect to the parameters.""" + target[0] += np.sum(partial) + + def dK_dX(self,X,X2,target): + """derivative of the covariance matrix with respect to X.""" + if X2 is None: X2 = X + dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscales),-1))[:,:,None] + ddist_dX = (X[:,None,:]-X2[None,:,:])/self.lengthscales**2/np.where(dist!=0.,dist,np.inf) + dK_dX += - np.transpose(3*self.variance*dist*np.exp(-np.sqrt(3)*dist)*ddist_dX,(1,0,2)) + target += np.sum(dK_dX*partial.T[:,:,None],0) + + def dKdiag_dX(self,X,target): + pass + def Gram_matrix(self,F,F1,F2,lower,upper): """ Return the Gram matrix of the vector of functions F with respect to the RKHS norm. The use of this function is limited to D=1. @@ -87,25 +116,3 @@ class Matern32(kernpart): #return(G) return(self.lengthscales**3/(12.*np.sqrt(3)*self.variance) * G + 1./self.variance*np.dot(Flower,Flower.T) + self.lengthscales**2/(3.*self.variance)*np.dot(F1lower,F1lower.T)) - def dK_dtheta(self,X,X2,target): - """derivative of the cross-covariance matrix with respect to the parameters (shape is NxMxNparam)""" - if X2 is None: X2 = X - dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscales),-1)) - dvar = (1+np.sqrt(3.)*dist)*np.exp(-np.sqrt(3.)*dist) - invdist = 1./np.where(dist!=0.,dist,np.inf) - dist2M = np.square(X[:,None,:]-X2[None,:,:])/self.lengthscales**3 - dl = (self.variance* 3 * dist * np.exp(-np.sqrt(3.)*dist))[:,:,np.newaxis] * dist2M*invdist[:,:,np.newaxis] - np.add(target[:,:,0],dvar, target[:,:,0]) - np.add(target[:,:,1:],dl, target[:,:,1:]) - def dKdiag_dtheta(self,X,target): - """derivative of the diagonal of the covariance matrix with respect to the parameters (shape is NxNparam)""" - np.add(target[:,0],1.,target[:,0]) - def dK_dX(self,X,X2,target): - """derivative of the covariance matrix with respect to X (*! shape is NxMxD !*).""" - if X2 is None: X2 = X - dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscales),-1))[:,:,None] - ddist_dX = (X[:,None,:]-X2[None,:,:])/self.lengthscales**2/np.where(dist!=0.,dist,np.inf) - target += - np.transpose(3*self.variance*dist*np.exp(-np.sqrt(3)*dist)*ddist_dX,(1,0,2)) - def dKdiag_dX(self,X,target): - pass - diff --git a/GPy/kern/Matern52.py b/GPy/kern/Matern52.py index 4af65a89..65059a5b 100644 --- a/GPy/kern/Matern52.py +++ b/GPy/kern/Matern52.py @@ -33,33 +33,61 @@ class Matern52(kernpart): self.Nparam = self.D + 1 self.name = 'Mat52' self.set_param(np.hstack((variance,lengthscales))) - self._Z, self._mu, self._S = np.empty(shape=(3,1)) # cached versions of Z,mu,S + def get_param(self): """return the value of the parameters.""" return np.hstack((self.variance,self.lengthscales)) + def set_param(self,x): """set the value of the parameters.""" assert x.size==(self.D+1) self.variance = x[0] self.lengthscales = x[1:] - self.lengthscales2 = np.square(self.lengthscales) - self._Z, self._mu, self._S = np.empty(shape=(3,1)) # cached versions of Z,mu,S + def get_param_names(self): """return parameter names.""" if self.D==1: return ['variance','lengthscale'] else: return ['variance']+['lengthscale_%i'%i for i in range(self.lengthscales.size)] + def K(self,X,X2,target): """Compute the covariance matrix between X and X2.""" if X2 is None: X2 = X dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscales),-1)) np.add(self.variance*(1+np.sqrt(5.)*dist+5./3*dist**2)*np.exp(-np.sqrt(5.)*dist), target,target) + def Kdiag(self,X,target): """Compute the diagonal of the covariance matrix associated to X.""" np.add(target,self.variance,target) + def dK_dtheta(self,partial,X,X2,target): + """derivative of the covariance matrix with respect to the parameters.""" + if X2 is None: X2 = X + dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscales),-1)) + invdist = 1./np.where(dist!=0.,dist,np.inf) + dist2M = np.square(X[:,None,:]-X2[None,:,:])/self.lengthscales**3 + dvar = (1+np.sqrt(5.)*dist+5./3*dist**2)*np.exp(-np.sqrt(5.)*dist) + dl = (self.variance * 5./3 * dist * (1 + np.sqrt(5.)*dist ) * np.exp(-np.sqrt(5.)*dist))[:,:,np.newaxis] * dist2M*invdist[:,:,np.newaxis] + target[0] += np.sum(dvar*partial) + target[1:] += (dl*partial[:,:,None]).sum(0).sum(0) + + def dKdiag_dtheta(self,X,target): + """derivative of the diagonal of the covariance matrix with respect to the parameters.""" + target[0] += np.sum(partial) + + def dK_dX(self,X,X2,target): + """derivative of the covariance matrix with respect to X.""" + if X2 is None: X2 = X + dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscales),-1))[:,:,None] + ddist_dX = (X[:,None,:]-X2[None,:,:])/self.lengthscales**2/np.where(dist!=0.,dist,np.inf) + dK_dX += - np.transpose(self.variance*5./3*dist*(1+np.sqrt(5)*dist)*np.exp(-np.sqrt(5)*dist)*ddist_dX,(1,0,2)) + target += np.sum(dK_dX*partial.T[:,:,None],0) + + def dKdiag_dX(self,X,target): + pass + def Gram_matrix(self,F,F1,F2,F3,lower,upper): """ Return the Gram matrix of the vector of functions F with respect to the RKHS norm. The use of this function is limited to D=1. @@ -91,26 +119,5 @@ class Matern52(kernpart): orig2 = 3./5*self.lengthscales**2 * ( np.dot(F1lower,F1lower.T) + 1./8*np.dot(Flower,F2lower.T) + 1./8*np.dot(F2lower,Flower.T)) return(1./self.variance* (G_coef*G + orig + orig2)) - def dK_dtheta(self,X,X2,target): - """derivative of the cross-covariance matrix with respect to the parameters (shape is NxMxNparam)""" - if X2 is None: X2 = X - dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscales),-1)) - invdist = 1./np.where(dist!=0.,dist,np.inf) - dist2M = np.square(X[:,None,:]-X2[None,:,:])/self.lengthscales**3 - dvar = (1+np.sqrt(5.)*dist+5./3*dist**2)*np.exp(-np.sqrt(5.)*dist) - dl = (self.variance * 5./3 * dist * (1 + np.sqrt(5.)*dist ) * np.exp(-np.sqrt(5.)*dist))[:,:,np.newaxis] * dist2M*invdist[:,:,np.newaxis] - np.add(target[:,:,0],dvar, target[:,:,0]) - np.add(target[:,:,1:],dl, target[:,:,1:]) - def dKdiag_dtheta(self,X,target): - """derivative of the diagonal of the covariance matrix with respect to the parameters (shape is NxNparam)""" - np.add(target[:,0],1.,target[:,0]) - def dK_dX(self,X,X2,target): - """derivative of the covariance matrix with respect to X (*! shape is NxMxD !*).""" - if X2 is None: X2 = X - dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscales),-1))[:,:,None] - ddist_dX = (X[:,None,:]-X2[None,:,:])/self.lengthscales**2/np.where(dist!=0.,dist,np.inf) - target += - np.transpose(self.variance*5./3*dist*(1+np.sqrt(5)*dist)*np.exp(-np.sqrt(5)*dist)*ddist_dX,(1,0,2)) - def dKdiag_dX(self,X,target): - pass diff --git a/GPy/kern/exponential.py b/GPy/kern/exponential.py index 402ebd82..ba97881e 100644 --- a/GPy/kern/exponential.py +++ b/GPy/kern/exponential.py @@ -62,7 +62,7 @@ class exponential(kernpart): np.add(target,self.variance,target) def dK_dtheta(self,partial,X,X2,target): - """derivative of the cross-covariance matrix with respect to the parameters (shape is NxMxNparam)""" + """derivative of the covariance matrix with respect to the parameters.""" if X2 is None: X2 = X dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscales),-1)) invdist = 1./np.where(dist!=0.,dist,np.inf) @@ -73,12 +73,12 @@ class exponential(kernpart): target[1:] += (dl*partial[:,:,None]).sum(0).sum(0) def dKdiag_dtheta(self,partial,X,target): - """derivative of the diagonal of the covariance matrix with respect to the parameters (shape is NxNparam)""" + """derivative of the diagonal of the covariance matrix with respect to the parameters.""" #NB: derivative of diagonal elements wrt lengthscale is 0 target[0] += np.sum(partial) def dK_dX(self,X,X2,target): - """derivative of the covariance matrix with respect to X (*! shape is NxMxD !*).""" + """derivative of the covariance matrix with respect to X.""" if X2 is None: X2 = X dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscales),-1))[:,:,None] ddist_dX = (X[:,None,:]-X2[None,:,:])/self.lengthscales**2/np.where(dist!=0.,dist,np.inf) From 049de98d3b510357f6f391c23f8be34eea159c5c Mon Sep 17 00:00:00 2001 From: Nicolas Date: Mon, 3 Dec 2012 11:26:10 +0000 Subject: [PATCH 3/3] rbf_ARD now in the updated format for the computation of the derivatives (included for the psi-statistics, but not tested) --- GPy/kern/rbf_ARD.py | 63 +++++++++++++++++++++++++++------------------ 1 file changed, 38 insertions(+), 25 deletions(-) diff --git a/GPy/kern/rbf_ARD.py b/GPy/kern/rbf_ARD.py index f1e5f36a..1f90bb0a 100644 --- a/GPy/kern/rbf_ARD.py +++ b/GPy/kern/rbf_ARD.py @@ -30,6 +30,7 @@ class rbf_ARD(kernpart): def get_param(self): return np.hstack((self.variance,self.lengthscales)) + def set_param(self,x): assert x.size==(self.D+1) self.variance = x[0] @@ -37,61 +38,73 @@ class rbf_ARD(kernpart): self.lengthscales2 = np.square(self.lengthscales) #reset cached results self._Z, self._mu, self._S = np.empty(shape=(3,1)) # cached versions of Z,mu,S + def get_param_names(self): if self.D==1: return ['variance','lengthscale'] else: return ['variance']+['lengthscale_%i'%i for i in range(self.lengthscales.size)] + def K(self,X,X2,target): self._K_computations(X,X2) np.add(self.variance*self._K_dvar, target,target) + def Kdiag(self,X,target): np.add(target,self.variance,target) - def dK_dtheta(self,X,X2,target): - """Return shape is NxMx(Ntheta)""" + + def dK_dtheta(self,partial,X,X2,target): self._K_computations(X,X2) dl = self._K_dvar[:,:,None]*self.variance*self._K_dist2/self.lengthscales - np.add(target[:,:,0],self._K_dvar, target[:,:,0]) - np.add(target[:,:,1:],dl, target[:,:,1:]) + target[0] += np.sum(self._K_dvar*partial) + target[1:] += (dl*partial[:,:,None]).sum(0).sum(0) + def dKdiag_dtheta(self,X,target): - np.add(target[:,0],1.,target[:,0]) + target[0] += np.sum(partial) + def dK_dX(self,X,X2,target): self._K_computations(X,X2) dZ = self.variance*self._K_dvar[:,:,None]*self._K_dist/self.lengthscales2 - np.add(target,-dZ.transpose(1,0,2),target) + dK_dX = -dZ.transpose(1,0,2) + target += np.sum(dK_dX*partial.T[:,:,None],0) + + def dKdiag_dX(self,partial,X,target): + pass + def psi0(self,Z,mu,S,target): - np.add(target, self.variance, target) - def dpsi0_dtheta(self,Z,mu,S,target): - target[:,0] += 1. + target += self.variance + + def dpsi0_dtheta(self,partial,Z,mu,S,target): + target[0] += 1. + def dpsi0_dmuS(self,Z,mu,S,target_mu,target_S): pass + def psi1(self,Z,mu,S,target): - """Think N,M,Q """ self._psi_computations(Z,mu,S) np.add(target, self._psi1,target) - def dpsi1_dtheta(self,Z,mu,S,target): - """N,Q,(Ntheta)""" + def dpsi1_dtheta(self,partial,Z,mu,S,target): self._psi_computations(Z,mu,S) denom_deriv = S[:,None,:]/(self.lengthscales**3+self.lengthscales*S[:,None,:]) d_length = self._psi1[:,:,None]*(self.lengthscales*np.square(self._psi1_dist/(self.lengthscales2+S[:,None,:])) + denom_deriv) - target[:,:,0] += self._psi1/self.variance - target[:,:,1:] += d_length - def dpsi1_dZ(self,Z,mu,S,target): + target[0] += np.sum(partial*self._psi1/self.variance) + target[1:] += (d_length*partial[:,:,None]).sum(0).sum(0) + + def dpsi1_dZ(self,partial,Z,mu,S,target): self._psi_computations(Z,mu,S) np.add(target,-self._psi1[:,:,None]*self._psi1_dist/self.lengthscales2/self._psi1_denom,target) + target += np.sum(partial[:,:,None]*-self._psi1[:,:,None]*self._psi1_dist/self.lengthscales2/self._psi1_denom,0) - def dpsi1_dmuS(self,Z,mu,S,target_mu,target_S): + def dpsi1_dmuS(self,partial,Z,mu,S,target_mu,target_S): """return shapes are N,M,Q""" self._psi_computations(Z,mu,S) tmp = self._psi1[:,:,None]/self.lengthscales2/self._psi1_denom - np.add(target_mu,tmp*self._psi1_dist,target_mu) - np.add(target_S, 0.5*tmp*(self._psi1_dist_sq-1), target_S) + target_mu += np.sum(partial*tmp*self._psi1_dist,1) + target_S += np.sum(partial*0.5*tmp*(self._psi1_dist_sq-1),1) def psi2(self,Z,mu,S,target): - """shape N,M,M""" self._psi_computations(Z,mu,S) - np.add(target, self._psi2,target) + target += self._psi2.sum(0) #TODO: psi2 should be NxMxM (for het. noise) def dpsi2_dtheta(self,Z,mu,S,target): """Shape N,M,M,Ntheta""" @@ -99,21 +112,21 @@ class rbf_ARD(kernpart): d_var = np.sum(2.*self._psi2/self.variance,0) d_length = self._psi2[:,:,:,None]*(0.5*self._psi2_Zdist_sq*self._psi2_denom + 2.*self._psi2_mudist_sq + 2.*S[:,None,None,:]/self.lengthscales2)/(self.lengthscales*self._psi2_denom) d_length = d_length.sum(0) - target[:,:,0] += d_var - target[:,:,1:] += d_length + target[0] += np.sum(partial*d_var) + target[1:] += (d_length*partial[:,:,None]).sum(0).sum(0) def dpsi2_dZ(self,Z,mu,S,target): """Returns shape N,M,M,Q""" self._psi_computations(Z,mu,S) dZ = self._psi2[:,:,:,None]/self.lengthscales2*(-0.5*self._psi2_Zdist + self._psi2_mudist/self._psi2_denom) - target += dZ + target += np.sum(partial[None,:,:,None]*dZ,0).sum(1) def dpsi2_dmuS(self,Z,mu,S,target_mu,target_S): """Think N,M,M,Q """ self._psi_computations(Z,mu,S) tmp = self._psi2[:,:,:,None]/self.lengthscales2/self._psi2_denom - np.add(target_mu, -tmp*(2.*self._psi2_mudist),target_mu) #N,M,M,Q - np.add(target_S, tmp*(2.*self._psi2_mudist_sq-1), target_S) #N,M,M,Q + target_mu += (partial*-tmp*2.*self._psi2_mudist).sum(1).sum(1) + target_S += (partial*tmp*(2.*self._psi2_mudist_sq-1)).sum(1).sum(1) def _K_computations(self,X,X2): if not (np.all(X==self._X) and np.all(X2==self._X2)):