diff --git a/GPy/__init__.py b/GPy/__init__.py index 6993d5c2..381d6232 100644 --- a/GPy/__init__.py +++ b/GPy/__init__.py @@ -6,5 +6,6 @@ import kern import models import inference import util +import examples #import examples TODO: discuss! from core import priors diff --git a/GPy/core/model.py b/GPy/core/model.py index 063eaf7d..77e66600 100644 --- a/GPy/core/model.py +++ b/GPy/core/model.py @@ -80,19 +80,22 @@ class model(parameterised): for w in which: self.priors[w] = what - def get(self,name): + def get(self,name, return_names=False): """ - get a model parameter by name + Get a model parameter by name. The name is applied as a regular expression and all parameters that match that regular expression are returned. """ matches = self.grep_param_names(name) if len(matches): - return self._get_params()[matches] + if return_names: + return self._get_params()[matches], np.asarray(self._get_param_names())[matches].tolist() + else: + return self._get_params()[matches] else: raise AttributeError, "no parameter matches %s"%name def set(self,name,val): """ - Set a model parameter by name + Set model parameter(s) by name. The name is provided as a regular expression. All parameters matching that regular expression are set to ghe given value. """ matches = self.grep_param_names(name) if len(matches): @@ -102,6 +105,20 @@ class model(parameterised): else: raise AttributeError, "no parameter matches %s"%name + def get_gradient(self,name, return_names=False): + """ + Get model gradient(s) by name. The name is applied as a regular expression and all parameters that match that regular expression are returned. + """ + matches = self.grep_param_names(name) + if len(matches): + if return_names: + return self._log_likelihood_gradients()[matches], np.asarray(self._get_param_names())[matches].tolist() + else: + return self._log_likelihood_gradients()[matches] + else: + raise AttributeError, "no parameter matches %s"%name + + def log_prior(self): diff --git a/GPy/examples/regression.py b/GPy/examples/regression.py index 79763504..9444e899 100644 --- a/GPy/examples/regression.py +++ b/GPy/examples/regression.py @@ -17,10 +17,8 @@ def toy_rbf_1d(): # create simple GP model m = GPy.models.GP_regression(data['X'],data['Y']) - # contrain all parameters to be positive - m.constrain_positive('') - # optimize + m.ensure_default_constraints() m.optimize() # plot @@ -35,10 +33,8 @@ def rogers_girolami_olympics(): # create simple GP model m = GPy.models.GP_regression(data['X'],data['Y']) - # contrain all parameters to be positive - m.constrain_positive('') - # optimize + m.ensure_default_constraints() m.optimize() # plot @@ -57,10 +53,8 @@ def toy_rbf_1d_50(): # create simple GP model m = GPy.models.GP_regression(data['X'],data['Y']) - # contrain all parameters to be positive - m.constrain_positive('') - # optimize + m.ensure_default_constraints() m.optimize() # plot @@ -75,10 +69,8 @@ def silhouette(): # create simple GP model m = GPy.models.GP_regression(data['X'],data['Y']) - # contrain all parameters to be positive - m.constrain_positive('') - # optimize + m.ensure_default_constraints() m.optimize() print(m) @@ -118,20 +110,15 @@ def multiple_optima(gene_number=937,resolution=80, model_restarts=10, seed=10000 kern = GPy.kern.rbf(1, variance=np.random.exponential(1.), lengthscale=np.random.exponential(50.)) + GPy.kern.white(1,variance=np.random.exponential(1.)) m = GPy.models.GP_regression(data['X'],data['Y'], kernel=kern) - params = m._get_params() - optim_point_x[0] = params[1] - optim_point_y[0] = np.log10(params[0]) - np.log10(params[2]); - - # contrain all parameters to be positive - m.constrain_positive('') + optim_point_x[0] = m.get('rbf_lengthscale') + optim_point_y[0] = np.log10(m.get('rbf_variance')) - np.log10(m.get('white_variance')); # optimize + m.ensure_default_constraints() m.optimize(xtol=1e-6,ftol=1e-6) - params = m._get_params() - optim_point_x[1] = params[1] - optim_point_y[1] = np.log10(params[0]) - np.log10(params[2]); - print(m) + optim_point_x[1] = m.get('rbf_lengthscale') + optim_point_y[1] = np.log10(m.get('rbf_variance')) - np.log10(m.get('white_variance')); pb.arrow(optim_point_x[0], optim_point_y[0], optim_point_x[1]-optim_point_x[0], optim_point_y[1]-optim_point_y[0], label=str(i), head_length=1, head_width=0.5, fc='k', ec='k') models.append(m) diff --git a/GPy/examples/sparse_GPLVM_demo.py b/GPy/examples/sparse_GPLVM_demo.py index 6ca6c941..3f1969fe 100644 --- a/GPy/examples/sparse_GPLVM_demo.py +++ b/GPy/examples/sparse_GPLVM_demo.py @@ -10,11 +10,11 @@ print "sparse GPLVM with RBF kernel" N = 100 M = 4 -Q = 1 +Q = 2 D = 2 #generate GPLVM-like data X = np.random.rand(N, Q) -k = GPy.kern.rbf(Q, 1.0, 2.0) + GPy.kern.white(Q, 0.00001) +k = GPy.kern.rbf(Q,1.,2*np.ones((1,))) + GPy.kern.white(Q, 0.00001) K = k.K(X) Y = np.random.multivariate_normal(np.zeros(N),K,D).T diff --git a/GPy/kern/Matern32.py b/GPy/kern/Matern32.py index 6517ac2c..1270e3f9 100644 --- a/GPy/kern/Matern32.py +++ b/GPy/kern/Matern32.py @@ -20,43 +20,52 @@ class Matern32(kernpart): :type D: int :param variance: the variance :math:`\sigma^2` :type variance: float - :param lengthscale: the lengthscales :math:`\ell_i` + :param lengthscale: the lengthscale :math:`\ell_i` :type lengthscale: np.ndarray of size (D,) :rtype: kernel object """ - def __init__(self,D,variance=1.,lengthscales=None): + def __init__(self,D,variance=1.,lengthscale=None,ARD=False): self.D = D - if lengthscales is not None: - assert lengthscales.shape==(self.D,) + self.ARD = ARD + if ARD == False: + self.Nparam = 2 + self.name = 'Mat32' + if lengthscale is not None: + assert lengthscale.shape == (1,) + else: + lengthscale = np.ones(1) else: - lengthscales = np.ones(self.D) - self.Nparam = self.D + 1 - self.name = 'Mat32' - self._set_params(np.hstack((variance,lengthscales))) + self.Nparam = self.D + 1 + self.name = 'Mat32_ARD' + if lengthscale is not None: + assert lengthscale.shape == (self.D,) + else: + lengthscale = np.ones(self.D) + self._set_params(np.hstack((variance,lengthscale))) def _get_params(self): """return the value of the parameters.""" - return np.hstack((self.variance,self.lengthscales)) + return np.hstack((self.variance,self.lengthscale)) def _set_params(self,x): """set the value of the parameters.""" - assert x.size==(self.D+1) + assert x.size == self.Nparam self.variance = x[0] - self.lengthscales = x[1:] + self.lengthscale = x[1:] def _get_param_names(self): """return parameter names.""" - if self.D==1: + if self.Nparam == 2: return ['variance','lengthscale'] else: - return ['variance']+['lengthscale_%i'%i for i in range(self.lengthscales.size)] + return ['variance']+['lengthscale_%i'%i for i in range(self.lengthscale.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)) + dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscale),-1)) np.add(self.variance*(1+np.sqrt(3.)*dist)*np.exp(-np.sqrt(3.)*dist), target,target) def Kdiag(self,X,target): @@ -66,13 +75,20 @@ class Matern32(kernpart): 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)) + dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscale),-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] + dist2M = np.square(X[:,None,:]-X2[None,:,:])/self.lengthscale**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) + if self.ARD == True: + dl = (self.variance* 3 * dist * np.exp(-np.sqrt(3.)*dist))[:,:,np.newaxis] * dist2M*invdist[:,:,np.newaxis] + #dl = self.variance*dvar[:,:,None]*dist2M*invdist[:,:,None] + target[1:] += (dl*partial[:,:,None]).sum(0).sum(0) + else: + dl = (self.variance* 3 * dist * np.exp(-np.sqrt(3.)*dist)) * dist2M.sum(-1)*invdist + #dl = self.variance*dvar*dist2M.sum(-1)*invdist + target[1] += np.sum(dl*partial) def dKdiag_dtheta(self,partial,X,target): """derivative of the diagonal of the covariance matrix with respect to the parameters.""" @@ -81,8 +97,8 @@ class Matern32(kernpart): def dK_dX(self,partial,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) + dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscale),-1))[:,:,None] + ddist_dX = (X[:,None,:]-X2[None,:,:])/self.lengthscale**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) @@ -104,7 +120,7 @@ class Matern32(kernpart): """ assert self.D == 1 def L(x,i): - return(3./self.lengthscales**2*F[i](x) + 2*np.sqrt(3)/self.lengthscales*F1[i](x) + F2[i](x)) + return(3./self.lengthscale**2*F[i](x) + 2*np.sqrt(3)/self.lengthscale*F1[i](x) + F2[i](x)) n = F.shape[0] G = np.zeros((n,n)) for i in range(n): @@ -114,5 +130,5 @@ class Matern32(kernpart): F1lower = np.array([f(lower) for f in F1])[:,None] #print "OLD \n", np.dot(F1lower,F1lower.T), "\n \n" #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)) + 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)) diff --git a/GPy/kern/__init__.py b/GPy/kern/__init__.py index cd893bac..4a36d6d0 100644 --- a/GPy/kern/__init__.py +++ b/GPy/kern/__init__.py @@ -2,5 +2,5 @@ # Licensed under the BSD 3-clause license (see LICENSE.txt) -from constructors import rbf, Matern32, Matern52, exponential, linear, white, bias, finite_dimensional, rbf_ARD, spline, Brownian, linear_ARD, rbf_sympy, sympykern +from constructors import rbf, Matern32, Matern52, exponential, linear, white, bias, finite_dimensional, spline, Brownian, linear_ARD, rbf_sympy, sympykern from kern import kern diff --git a/GPy/kern/constructors.py b/GPy/kern/constructors.py index 0ddc09e3..5f676d9b 100644 --- a/GPy/kern/constructors.py +++ b/GPy/kern/constructors.py @@ -22,7 +22,7 @@ from Brownian import Brownian as Brownianpart #using meta-classes to make the objects construct properly wthout them. -def rbf(D,variance=1., lengthscale=1.): +def rbf(D,variance=1., lengthscale=None,ARD=False): """ Construct an RBF kernel @@ -33,21 +33,7 @@ def rbf(D,variance=1., lengthscale=1.): :param lengthscale: the lengthscale of the kernel :type lengthscale: float """ - part = rbfpart(D,variance,lengthscale) - return kern(D, [part]) - -def rbf_ARD(D,variance=1., lengthscales=None): - """ - Construct an RBF kernel with Automatic Relevance Determination (ARD) - - :param D: dimensionality of the kernel, obligatory - :type D: int - :param variance: the variance of the kernel - :type variance: float - :param lengthscales: the lengthscales of the kernel - :type lengthscales: None|np.ndarray - """ - part = rbf_ARD_part(D,variance,lengthscales) + part = rbfpart(D,variance,lengthscale,ARD) return kern(D, [part]) def linear(D,lengthscales=None): @@ -86,7 +72,7 @@ def white(D,variance=1.): part = whitepart(D,variance) return kern(D, [part]) -def exponential(D,variance=1., lengthscales=None): +def exponential(D,variance=1., lengthscale=None, ARD=False): """ Construct a exponential kernel. @@ -96,10 +82,10 @@ def exponential(D,variance=1., lengthscales=None): variance (float) lengthscales (np.ndarray) """ - part = exponentialpart(D,variance, lengthscales) + part = exponentialpart(D,variance, lengthscale, ARD) return kern(D, [part]) -def Matern32(D,variance=1., lengthscales=None): +def Matern32(D,variance=1., lengthscale=None, ARD=False): """ Construct a Matern 3/2 kernel. @@ -109,7 +95,7 @@ def Matern32(D,variance=1., lengthscales=None): variance (float) lengthscales (np.ndarray) """ - part = Matern32part(D,variance, lengthscales) + part = Matern32part(D,variance, lengthscale, ARD) return kern(D, [part]) def Matern52(D,variance=1., lengthscales=None): diff --git a/GPy/kern/exponential.py b/GPy/kern/exponential.py index 2df6a958..0ea1e922 100644 --- a/GPy/kern/exponential.py +++ b/GPy/kern/exponential.py @@ -24,37 +24,46 @@ class exponential(kernpart): :rtype: kernel object """ - def __init__(self,D,variance=1.,lengthscales=None): + def __init__(self,D,variance=1.,lengthscale=None,ARD=False): self.D = D - if lengthscales is not None: - assert lengthscales.shape==(self.D,) + self.ARD = ARD + if ARD == False: + self.Nparam = 2 + self.name = 'exp' + if lengthscale is not None: + assert lengthscale.shape == (1,) + else: + lengthscale = np.ones(1) else: - lengthscales = np.ones(self.D) - self.Nparam = self.D + 1 - self.name = 'exp' - self._set_params(np.hstack((variance,lengthscales))) + self.Nparam = self.D + 1 + self.name = 'exp_ARD' + if lengthscale is not None: + assert lengthscale.shape == (self.D,) + else: + lengthscale = np.ones(self.D) + self._set_params(np.hstack((variance,lengthscale))) def _get_params(self): """return the value of the parameters.""" - return np.hstack((self.variance,self.lengthscales)) + return np.hstack((self.variance,self.lengthscale)) def _set_params(self,x): """set the value of the parameters.""" - assert x.size==(self.D+1) + assert x.size == self.Nparam self.variance = x[0] - self.lengthscales = x[1:] + self.lengthscale = x[1:] def _get_param_names(self): """return parameter names.""" - if self.D==1: + if self.Nparam == 2: return ['variance','lengthscale'] else: - return ['variance']+['lengthscale_%i'%i for i in range(self.lengthscales.size)] + return ['variance']+['lengthscale_%i'%i for i in range(self.lengthscale.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)) + dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscale),-1)) np.add(self.variance*np.exp(-dist), target,target) def Kdiag(self,X,target): @@ -64,13 +73,17 @@ class exponential(kernpart): 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)) + dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscale),-1)) invdist = 1./np.where(dist!=0.,dist,np.inf) - dist2M = np.square(X[:,None,:]-X2[None,:,:])/self.lengthscales**3 + dist2M = np.square(X[:,None,:]-X2[None,:,:])/self.lengthscale**3 dvar = np.exp(-dist) - dl = self.variance*dvar[:,:,None]*dist2M*invdist[:,:,None] target[0] += np.sum(dvar*partial) - target[1:] += (dl*partial[:,:,None]).sum(0).sum(0) + if self.ARD == True: + dl = self.variance*dvar[:,:,None]*dist2M*invdist[:,:,None] + target[1:] += (dl*partial[:,:,None]).sum(0).sum(0) + else: + dl = self.variance*dvar*dist2M.sum(-1)*invdist + target[1] += np.sum(dl*partial) def dKdiag_dtheta(self,partial,X,target): """derivative of the diagonal of the covariance matrix with respect to the parameters.""" @@ -80,8 +93,8 @@ class exponential(kernpart): def dK_dX(self,partial,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) + dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscale),-1))[:,:,None] + ddist_dX = (X[:,None,:]-X2[None,:,:])/self.lengthscale**2/np.where(dist!=0.,dist,np.inf) dK_dX = - np.transpose(self.variance*np.exp(-dist)*ddist_dX,(1,0,2)) target += np.sum(dK_dX*partial.T[:,:,None],0) @@ -101,14 +114,14 @@ class exponential(kernpart): """ assert self.D == 1 def L(x,i): - return(1./self.lengthscales*F[i](x) + F1[i](x)) + return(1./self.lengthscale*F[i](x) + F1[i](x)) n = F.shape[0] G = np.zeros((n,n)) for i in range(n): for j in range(i,n): G[i,j] = G[j,i] = integrate.quad(lambda x : L(x,i)*L(x,j),lower,upper)[0] Flower = np.array([f(lower) for f in F])[:,None] - return(self.lengthscales/2./self.variance * G + 1./self.variance * np.dot(Flower,Flower.T)) + return(self.lengthscale/2./self.variance * G + 1./self.variance * np.dot(Flower,Flower.T)) diff --git a/GPy/kern/rbf.py b/GPy/kern/rbf.py index c929369f..13ecaf86 100644 --- a/GPy/kern/rbf.py +++ b/GPy/kern/rbf.py @@ -20,17 +20,33 @@ class rbf(kernpart): :type D: int :param variance: the variance of the kernel :type variance: float - :param lengthscale: the lengthscale of the kernel - :type lengthscale: float + :param lengthscale: the vector of lengthscale of the kernel + :type lengthscale: np.ndarray + :param ARD: Auto Relevance Determination. If equal to "False", the kernel is isotropic (ie. one single lengthscale parameter \ell), otherwise there is one lengthscale parameter per dimension. + :type ARD: Boolean - .. Note: for rbf with different lengthscale on each dimension, see rbf_ARD """ - def __init__(self,D,variance=1.,lengthscale=1.): + def __init__(self,D,variance=1.,lengthscale=None,ARD=False): self.D = D - self.Nparam = 2 - self.name = 'rbf' - self._set_params(np.hstack((variance,lengthscale))) + self.ARD = ARD + if ARD == False: + self.Nparam = 2 + self.name = 'rbf' + if lengthscale is not None: + assert lengthscale.shape == (1,) + else: + lengthscale = np.ones(1) + + else: + self.Nparam = self.D + 1 + self.name = 'rbf_ARD' + if lengthscale is not None: + assert lengthscale.shape == (self.D,) + else: + lengthscale = np.ones(self.D) + + self._set_params(np.hstack((variance,lengthscale))) #initialize cache self._Z, self._mu, self._S = np.empty(shape=(3,1)) @@ -40,14 +56,19 @@ class rbf(kernpart): return np.hstack((self.variance,self.lengthscale)) def _set_params(self,x): - self.variance, self.lengthscale = x + assert x.size==(self.Nparam) + self.variance = x[0] + self.lengthscale = x[1:] self.lengthscale2 = np.square(self.lengthscale) #reset cached results self._X, self._X2, self._params = np.empty(shape=(3,1)) self._Z, self._mu, self._S = np.empty(shape=(3,1)) # cached versions of Z,mu,S def _get_param_names(self): - return ['variance','lengthscale'] + if self.Nparam == 2: + return ['variance','lengthscale'] + else: + return ['variance']+['lengthscale_%i'%i for i in range(self.lengthscale.size)] def K(self,X,X2,target): if X2 is None: @@ -61,7 +82,12 @@ class rbf(kernpart): def dK_dtheta(self,partial,X,X2,target): self._K_computations(X,X2) target[0] += np.sum(self._K_dvar*partial) - target[1] += np.sum(self._K_dvar*self.variance*self._K_dist2/self.lengthscale*partial) + if self.ARD == True: + dl = self._K_dvar[:,:,None]*self.variance*self._K_dist2/self.lengthscale + target[1:] += (dl*partial[:,:,None]).sum(0).sum(0) + else: + target[1] += np.sum(self._K_dvar*self.variance*(self._K_dist2.sum(-1))/self.lengthscale*partial) + #np.sum(self._K_dvar*self.variance*self._K_dist2/self.lengthscale*partial) def dKdiag_dtheta(self,partial,X,target): #NB: derivative of diagonal elements wrt lengthscale is 0 @@ -81,15 +107,13 @@ class rbf(kernpart): self._X = X self._X2 = X2 if X2 is None: X2 = X - XXT = np.dot(X,X2.T) - if X is X2: - self._K_dist2 = (-2.*XXT + np.diag(XXT)[:,np.newaxis] + np.diag(XXT)[np.newaxis,:])/self.lengthscale2 - else: - self._K_dist2 = (-2.*XXT + np.sum(np.square(X),1)[:,None] + np.sum(np.square(X2),1)[None,:])/self.lengthscale2 - # TODO Remove comments if this is fine. - # Commented out by Neil as doesn't seem to be used elsewhere. - #self._K_exponent = -0.5*self._K_dist2 - self._K_dvar = np.exp(-0.5*self._K_dist2) + self._K_dist = X[:,None,:]-X2[None,:,:] # this can be computationally heavy + self._params = np.empty(shape=(1,0))#ensure the next section gets called + if not np.all(self._params == self._get_params()): + self._params == self._get_params() + self._K_dist2 = np.square(self._K_dist/self.lengthscale) + #self._K_exponent = -0.5*self._K_dist2.sum(-1) #ND: commented out because seems not to be used + self._K_dvar = np.exp(-0.5*self._K_dist2.sum(-1)) def psi0(self,Z,mu,S,target): target += self.variance @@ -132,7 +156,7 @@ class rbf(kernpart): d_length = self._psi2[:,:,:,None]*(0.5*self._psi2_Zdist_sq*self._psi2_denom + 2.*self._psi2_mudist_sq + 2.*S[:,None,None,:]/self.lengthscale2)/(self.lengthscale*self._psi2_denom) d_length = d_length.sum(0) target[0] += np.sum(partial*d_var) - target[1] += np.sum(d_length*partial) + target[1:] += (d_length*partial[:,:,None]).sum(0).sum(0) def dpsi2_dZ(self,partial,Z,mu,S,target): """Returns shape N,M,M,Q""" @@ -175,4 +199,3 @@ class rbf(kernpart): self._psi2 = np.square(self.variance)*np.exp(self._psi2_exponent) # N,M,M self._Z, self._mu, self._S = Z, mu,S - diff --git a/GPy/models/GP_regression.py b/GPy/models/GP_regression.py index eee0fb58..72a24307 100644 --- a/GPy/models/GP_regression.py +++ b/GPy/models/GP_regression.py @@ -63,10 +63,10 @@ class GP_regression(model): self._Ystd = np.ones((1,self.Y.shape[1])) if self.D > self.N: - # then it's more efficient to store Youter - self.Youter = np.dot(self.Y, self.Y.T) + # then it's more efficient to store YYT + self.YYT = np.dot(self.Y, self.Y.T) else: - self.Youter = None + self.YYT = None model.__init__(self) @@ -83,23 +83,23 @@ class GP_regression(model): def _model_fit_term(self): """ - Computes the model fit using Youter if it's available + Computes the model fit using YYT if it's available """ - if self.Youter is None: + if self.YYT is None: return -0.5*np.sum(np.square(np.dot(self.Li,self.Y))) else: - return -0.5*np.sum(np.multiply(self.Ki, self.Youter)) + return -0.5*np.sum(np.multiply(self.Ki, self.YYT)) def log_likelihood(self): complexity_term = -0.5*self.N*self.D*np.log(2.*np.pi) - 0.5*self.D*self.K_logdet return complexity_term + self._model_fit_term() def dL_dK(self): - if self.Youter is None: + if self.YYT is None: alpha = np.dot(self.Ki,self.Y) dL_dK = 0.5*(np.dot(alpha,alpha.T)-self.D*self.Ki) else: - dL_dK = 0.5*(mdot(self.Ki, self.Youter, self.Ki) - self.D*self.Ki) + dL_dK = 0.5*(mdot(self.Ki, self.YYT, self.Ki) - self.D*self.Ki) return dL_dK diff --git a/GPy/models/generalized_FITC.py b/GPy/models/generalized_FITC.py index d8e9c23d..a5ed8d0a 100644 --- a/GPy/models/generalized_FITC.py +++ b/GPy/models/generalized_FITC.py @@ -91,9 +91,9 @@ class generalized_FITC(model): def log_likelihood(self): self.posterior_param() - self.Youter = np.dot(self.mu_tilde,self.mu_tilde.T) + self.YYT = np.dot(self.mu_tilde,self.mu_tilde.T) A = -self.hld - B = -.5*np.sum(self.Qi*self.Youter) + B = -.5*np.sum(self.Qi*self.YYT) C = sum(np.log(self.ep_approx.Z_hat)) D = .5*np.sum(np.log(1./self.ep_approx.tau_tilde + 1./self.ep_approx.tau_)) E = .5*np.sum((self.ep_approx.v_/self.ep_approx.tau_ - self.mu_tilde.flatten())**2/(1./self.ep_approx.tau_ + 1./self.ep_approx.tau_tilde)) diff --git a/GPy/models/warped_GP.py b/GPy/models/warped_GP.py index ff6267d2..8ce80c76 100644 --- a/GPy/models/warped_GP.py +++ b/GPy/models/warped_GP.py @@ -48,9 +48,9 @@ class warpedGP(GP_regression): # this supports the 'smart' behaviour in GP_regression if self.D > self.N: - self.Youter = np.dot(self.Y, self.Y.T) + self.YYT = np.dot(self.Y, self.Y.T) else: - self.Youter = None + self.YYT = None return self.Y diff --git a/GPy/testing/unit_tests.py b/GPy/testing/unit_tests.py index ff9aba0e..d2ef87f7 100644 --- a/GPy/testing/unit_tests.py +++ b/GPy/testing/unit_tests.py @@ -121,7 +121,7 @@ class GradientTests(unittest.TestCase): """ Testing GPLVM with rbf + bias and white kernel """ N, Q, D = 50, 1, 2 X = np.random.rand(N, Q) - k = GPy.kern.rbf(Q, 0.5, 0.9) + GPy.kern.bias(Q, 0.1) + GPy.kern.white(Q, 0.05) + k = GPy.kern.rbf(Q, 0.5, 0.9*np.ones((1,))) + GPy.kern.bias(Q, 0.1) + GPy.kern.white(Q, 0.05) K = k.K(X) Y = np.random.multivariate_normal(np.zeros(N),K,D).T m = GPy.models.GPLVM(Y, Q, kernel = k)