From d8eb155622e51a4a4fb62118af696b7d57c21aa8 Mon Sep 17 00:00:00 2001 From: Ricardo Andrade Date: Wed, 30 Jan 2013 16:00:03 +0000 Subject: [PATCH] Working for regression, still some bugs for EP. --- GPy/examples/sparse_ep_fix.py | 31 +++++++++++++------- GPy/models/sparse_GP.py | 55 ++++++++++++++++++----------------- 2 files changed, 50 insertions(+), 36 deletions(-) diff --git a/GPy/examples/sparse_ep_fix.py b/GPy/examples/sparse_ep_fix.py index ff90f2bb..defcb4eb 100644 --- a/GPy/examples/sparse_ep_fix.py +++ b/GPy/examples/sparse_ep_fix.py @@ -32,18 +32,29 @@ noise = GPy.kern.white(1) kernel = rbf + noise # create simple GP model -m = GPy.models.sparse_GP(X,Y=None, kernel=kernel, M=M,likelihood= likelihood) -#m = GPy.models.sparse_GP(X, Y, kernel, M=M) +#m = GPy.models.sparse_GP(X,Y=None, kernel=kernel, M=M,likelihood= likelihood) # contrain all parameters to be positive +#m.constrain_fixed('prec',100.) +m = GPy.models.sparse_GP(X, Y, kernel, M=M) m.ensure_default_constraints() -if not isinstance(m.likelihood,GPy.inference.likelihoods.gaussian): - m.approximate_likelihood() +#if not isinstance(m.likelihood,GPy.inference.likelihoods.gaussian): +# m.approximate_likelihood() print m.checkgrad() -#check gradient FIXME unit test please -# optimize and plot -#m.optimize('tnc', messages = 1) -m.EM() -m.plot(samples=3,full_cov=False) -# print(m) +m.optimize('tnc', messages = 1) +m.plot(samples=3) +print m +n = GPy.models.sparse_GP(X,Y=None, kernel=kernel, M=M,likelihood= likelihood) +n.ensure_default_constraints() +if not isinstance(n.likelihood,GPy.inference.likelihoods.gaussian): + n.approximate_likelihood() +print n.checkgrad() +pb.figure() +n.plot() + +""" +m = GPy.models.sparse_GP_regression(X, Y, kernel, M=M) +m.ensure_default_constraints() +print m.checkgrad() +""" diff --git a/GPy/models/sparse_GP.py b/GPy/models/sparse_GP.py index ba07254f..7f287174 100644 --- a/GPy/models/sparse_GP.py +++ b/GPy/models/sparse_GP.py @@ -10,6 +10,7 @@ from GP import GP from ..inference.EP import Full,DTC,FITC from ..inference.likelihoods import likelihood,probit,poisson,gaussian + #Still TODO: # make use of slices properly (kernel can now do this) # enable heteroscedatic noise (kernel will need to compute psi2 as a (NxMxM) array) @@ -35,12 +36,6 @@ class sparse_GP(GP): :type beta: float :param normalize_(X|Y) : whether to normalize the data before computing (predictions will be in original scales) :type normalize_(X|Y): bool - :parm likelihood: a GPy likelihood, defaults to gaussian - :param method_ep: sparse approximation used by Expectation Propagation algorithm, defaults to DTC - :type M: string (Full|DTC|FITC) - :param epsilon_ep: convergence criterion for the Expectation Propagation algorithm, defaults to 0.1 - :param powerep: power-EP parameters [$\eta$,$\delta$], defaults to [1.,1.] - :type powerep: list """ def __init__(self,X,Y=None,kernel=None,X_uncertainty=None,beta=100.,Z=None,Zslices=None,M=10,normalize_X=False,normalize_Y=False,likelihood=None,method_ep='DTC',epsilon_ep=1e-3,power_ep=[1.,1.]): @@ -70,20 +65,21 @@ class sparse_GP(GP): else: self.method_ep = method_ep + #normalise X uncertainty also + if self.has_uncertain_inputs: + self.X_uncertainty /= np.square(self._Xstd) def _set_params(self, p): self.Z = p[:self.M*self.Q].reshape(self.M, self.Q) if not self.EP: self.beta = p[self.M*self.Q] - #self.beta = np.repeat(p[self.M*self.Q],self.N)[:,None] self.kern._set_params(p[self.Z.size + 1:]) - self.beta2 = self.beta**2 else: self.kern._set_params(p[self.Z.size:]) if self.Y is None: self.Y = np.ones([self.N,1]) self._compute_kernel_matrices() - self._computations() #NOTE At this point computations of dL are not needed + self._computations() def _get_params(self): if not self.EP: @@ -97,19 +93,22 @@ class sparse_GP(GP): else: return sum([['iip_%i_%i'%(i,j) for i in range(self.Z.shape[0])] for j in range(self.Z.shape[1])],[]) + self.kern._get_param_names_transformed() + def _compute_kernel_matrices(self): # kernel computations, using BGPLVM notation #TODO: slices for psi statistics (easy enough) + self.Kmm = self.kern.K(self.Z) if self.has_uncertain_inputs: if not self.EP: - self.psi0 = self.kern.psi0(self.Z,self.X, self.X_uncertainty)#.sum() + self.psi0 = self.kern.psi0(self.Z,self.X, self.X_uncertainty)#.sum() NOTE psi0 is now a vector self.psi1 = self.kern.psi1(self.Z,self.X, self.X_uncertainty).T - self.psi2 = self.kern.psi2(self.Z,self.X, self.X_uncertainty)#FIXME add beta vector + self.psi2 = self.kern.psi2(self.Z,self.X, self.X_uncertainty) + #self.psi2_beta_scaled = ? else: raise NotImplementedError, "uncertain_inputs not yet supported for EP" else: - self.psi0 = self.kern.Kdiag(self.X,slices=self.Xslices)#.sum() FIXME + self.psi0 = self.kern.Kdiag(self.X,slices=self.Xslices)#.sum() self.psi1 = self.kern.K(self.Z,self.X) self.psi2 = np.dot(self.psi1,self.psi1.T) self.psi2_beta_scaled = np.dot(self.psi1,self.beta*self.psi1.T) @@ -124,22 +123,29 @@ class sparse_GP(GP): self.B = np.eye(self.M) + self.A self.Bi, self.LB, self.LBi, self.B_logdet = pdinv(self.B) self.LLambdai = np.dot(self.LBi, self.Lmi) - self.trace_K = self.psi0.sum() - np.trace(self.A) self.LBL_inv = mdot(self.Lmi.T, self.Bi, self.Lmi) self.C = mdot(self.LLambdai, self.psi1V) self.G = mdot(self.LBL_inv, self.psi1VVpsi1, self.LBL_inv.T) + self.trace_K_beta_scaled = (self.psi0*self.beta).sum() - np.trace(self.A) + if not self.EP: + self.trace_K = self.psi0.sum() - np.trace(self.A)/self.beta # Compute dL_dpsi - self.dL_dpsi0 = - 0.5 * self.D * self.beta.flatten() * np.ones(self.N) self.dL_dpsi1 = mdot(self.LLambdai.T,self.C,self.V.T) - #self.dL_dpsi2 = - 0.5 * self.beta * (self.D*(self.LBL_inv - self.Kmmi) + self.G) - self.dL_dpsi2 = - 0.5 * (self.D*(self.LBL_inv - self.Kmmi) + self.G) + if not self.EP: + self.dL_dpsi0 = - 0.5 * self.D * self.beta * np.ones(self.N) + if self.has_uncertain_inputs: + self.dL_dpsi2 = - 0.5 * self.beta * (self.D*(self.LBL_inv - self.Kmmi) + self.G) + else: + self.dL_dpsi2_ = - 0.5 * (self.D*(self.LBL_inv - self.Kmmi) + self.G) + else: + self.dL_dpsi0 = - 0.5 * self.D * self.beta.flatten() + if not self.has_uncertain_inputs: + self.dL_dpsi2_ = - 0.5 * (self.D*(self.LBL_inv - self.Kmmi) + self.G) # Compute dL_dKmm self.dL_dKmm = -0.5 * self.D * mdot(self.Lmi.T, self.A, self.Lmi) # dB - #self.dL_dKmm += -0.5 * self.D * (- self.LBL_inv - 2.*self.beta*mdot(self.LBL_inv, self.psi2, self.Kmmi) + self.Kmmi) # dC self.dL_dKmm += -0.5 * self.D * (- self.LBL_inv - 2.*mdot(self.LBL_inv, self.psi2_beta_scaled, self.Kmmi) + self.Kmmi) # dC - #self.dL_dKmm += np.dot(np.dot(self.G,self.beta*self.psi2) - np.dot(self.LBL_inv, self.psi1VVpsi1), self.Kmmi) + 0.5*self.G # dE self.dL_dKmm += np.dot(np.dot(self.G,self.psi2_beta_scaled) - np.dot(self.LBL_inv, self.psi1VVpsi1), self.Kmmi) + 0.5*self.G # dE def approximate_likelihood(self): @@ -164,7 +170,7 @@ class sparse_GP(GP): else: A = -0.5*self.D*(self.N*np.log(2.*np.pi) - np.sum(np.log(self.beta))) D = -0.5*self.trbetaYYT - B = -0.5*self.D*self.trace_K + B = -0.5*self.D*self.trace_K_beta_scaled C = -0.5*self.D * self.B_logdet E = +0.5*np.sum(self.psi1VVpsi1 * self.LBL_inv) return A+B+C+D+E @@ -194,7 +200,7 @@ class sparse_GP(GP): dL_dtheta += self.kern.dpsi2_dtheta(self.dL_dpsi2,self.Z,self.X, self.X_uncertainty) # for multiple_beta, dL_dpsi2 will be a different shape else: #re-cast computations in psi2 back to psi1: - dL_dpsi1 = self.dL_dpsi1 + 2.*np.dot(self.dL_dpsi2,self.psi1) + dL_dpsi1 = self.dL_dpsi1 + 2.*np.dot(self.dL_dpsi2_,self.beta.T*self.psi1) #dL_dpsi1 = self.dL_dpsi1 + 2.*np.dot(self.dL_dpsi2,self.psi1) dL_dtheta += self.kern.dK_dtheta(dL_dpsi1,self.Z,self.X) dL_dtheta += self.kern.dKdiag_dtheta(self.dL_dpsi0, self.X) @@ -210,7 +216,7 @@ class sparse_GP(GP): dL_dZ += self.kern.dpsi2_dZ(self.dL_dpsi2,self.Z,self.X, self.X_uncertainty) else: #re-cast computations in psi2 back to psi1: - dL_dpsi1 = self.dL_dpsi1 + 2.*np.dot(self.dL_dpsi2,self.psi1) + dL_dpsi1 = self.dL_dpsi1 + 2.*np.dot(self.dL_dpsi2_,self.beta.T*self.psi1)#dL_dpsi1 = self.dL_dpsi1 + 2.*np.dot(self.dL_dpsi2,self.psi1) dL_dZ += self.kern.dK_dX(dL_dpsi1,self.Z,self.X) return dL_dZ @@ -229,16 +235,14 @@ class sparse_GP(GP): Kxx = self.kern.K(Xnew) var = Kxx - mdot(Kx.T, (self.Kmmi - self.LBL_inv), Kx) if not self.EP: - var += np.eye(Xnew.shape[0])/self.beta # TODO: This beta doesn't belong here in the EP case. + var += np.eye(Xnew.shape[0])/self.beta else: raise NotImplementedError, "full_cov = True not implemented for EP" - #var = np.diag(var)[:,None] - #phi = self.likelihood.predictive_mean(mu,var) else: Kxx = self.kern.Kdiag(Xnew) var = Kxx - np.sum(Kx*np.dot(self.Kmmi - self.LBL_inv, Kx),0) if not self.EP: - var += 1./self.beta # TODO: This beta doesn't belong here in the EP case. + var += 1./self.beta else: phi = self.likelihood.predictive_mean(mu,var) return mu,var,phi @@ -247,7 +251,6 @@ class sparse_GP(GP): """ Plot the fitted model: just call the GP_regression plot function and then add inducing inputs """ - #GP_regression.plot(self,*args,**kwargs) GP.plot(self,*args,**kwargs) if self.Q==1: pb.plot(self.Z,self.Z*0+pb.ylim()[0],'k|',mew=1.5,markersize=12)