From 71b845eb603833eba01ea80d5eaa3a0493011c9a Mon Sep 17 00:00:00 2001 From: Ricardo Date: Wed, 8 May 2013 07:09:00 +0100 Subject: [PATCH] Some changes according to the changes in sparse_GP --- GPy/models/generalized_FITC.py | 49 ++++++++++++++++++++++++---------- 1 file changed, 35 insertions(+), 14 deletions(-) diff --git a/GPy/models/generalized_FITC.py b/GPy/models/generalized_FITC.py index 25b6c18f..966cbd39 100644 --- a/GPy/models/generalized_FITC.py +++ b/GPy/models/generalized_FITC.py @@ -9,6 +9,12 @@ from .. import kern from scipy import stats, linalg from sparse_GP import sparse_GP +def backsub_both_sides(L,X): + """ Return L^-T * X * L^-1, assumuing X is symmetrical and L is lower cholesky""" + tmp,_ = linalg.lapack.flapack.dtrtrs(L,np.asfortranarray(X),lower=1,trans=1) + return linalg.lapack.flapack.dtrtrs(L,np.asfortranarray(tmp.T),lower=1,trans=1)[0].T + + class generalized_FITC(sparse_GP): """ Naish-Guzman, A. and Holden, S. (2008) implemantation of EP with FITC. @@ -33,7 +39,7 @@ class generalized_FITC(sparse_GP): self.Z = Z self.M = self.Z.shape[0] - self._precision = likelihood.precision + self.true_precision = likelihood.precision sparse_GP.__init__(self, X, likelihood, kernel=kernel, Z=self.Z, X_variance=None, normalize_X=False) @@ -51,13 +57,16 @@ class generalized_FITC(sparse_GP): For a Gaussian (or direct: TODO) likelihood, no iteration is required: this function does nothing + + Diag(Knn - Qnn) is added to the noise term to use the tools already implemented in sparse_GP. + The true precison is now 'true_precision' not 'precision'. """ if self.has_uncertain_inputs: raise NotImplementedError, "FITC approximation not implemented for uncertain inputs" else: self.likelihood.fit_FITC(self.Kmm,self.psi1,self.psi0) - self._precision = self.likelihood.precision # Save the true precision - self.likelihood.precision = self._precision/(1. + self._precision*self.Diag0[:,None]) # Add the diagonal element of the FITC approximation + self.true_precision = self.likelihood.precision # Save the true precision + self.likelihood.precision = self.true_precision/(1. + self.true_precision*self.Diag0[:,None]) # Add the diagonal element of the FITC approximation self._set_params(self._get_params()) # update the GP def _FITC_computations(self): @@ -69,23 +78,23 @@ class generalized_FITC(sparse_GP): - removes the extra terms computed in the sparse_GP approximation - computes the likelihood gradients wrt the true precision. """ - #NOTE the true precison is now '_precison' not 'precision' + #NOTE the true precison is now 'true_precision' not 'precision' if self.likelihood.is_heteroscedastic: # Compute generalized FITC's diagonal term of the covariance - self.Qnn = mdot(self.psi1.T,self.Kmmi,self.psi1) + self.Lmi,info = linalg.lapack.flapack.dtrtrs(self.Lm,np.eye(self.M),lower=1) + Lmipsi1 = np.dot(self.Lmi,self.psi1) + self.Qnn = np.dot(Lmipsi1.T,Lmipsi1) + #self.Kmmi, Lm, Lmi, Kmm_logdet = pdinv(self.Kmm) + #self.Qnn = mdot(self.psi1.T,self.Kmmi,self.psi1) + #a = kj self.Diag0 = self.psi0 - np.diag(self.Qnn) - Iplus_Dprod_i = 1./(1.+ self.Diag0 * self._precision.flatten()) + Iplus_Dprod_i = 1./(1.+ self.Diag0 * self.true_precision.flatten()) self.Diag = self.Diag0 * Iplus_Dprod_i - #self.Diag = self.Diag0/(1.+ self.Diag0 * self._precision.flatten()) - self.P = Iplus_Dprod_i[:,None] * self.psi1.T - #self.P = (self.Diag / self.Diag0)[:,None] * self.psi1.T self.RPT0 = np.dot(self.Lmi,self.psi1) self.L = np.linalg.cholesky(np.eye(self.M) + np.dot(self.RPT0,((1. - Iplus_Dprod_i)/self.Diag0)[:,None]*self.RPT0.T)) - #self.L = np.linalg.cholesky(np.eye(self.M) + np.dot(self.RPT0,(1./self.Diag0 - Iplus_Dprod_i/self.Diag0)[:,None]*self.RPT0.T)) - #self.L = np.linalg.cholesky(np.eye(self.M) + np.dot(self.RPT0,(1./self.Diag0 - self.Diag/(self.Diag0**2))[:,None]*self.RPT0.T)) self.R,info = linalg.flapack.dtrtrs(self.L,self.Lmi,lower=1) self.RPT = np.dot(self.R,self.P.T) self.Sigma = np.diag(self.Diag) + np.dot(self.RPT.T,self.RPT) @@ -94,7 +103,16 @@ class generalized_FITC(sparse_GP): self.mu = self.w + np.dot(self.P,self.gamma) # Remove extra term from dL_dpsi1 - self.dL_dpsi1 -= mdot(self.Kmmi,self.psi1*self.likelihood.precision.flatten().reshape(1,self.N)) #dB + self.dL_dpsi1 -= mdot(self.Lmi.T,Lmipsi1*self.likelihood.precision.flatten().reshape(1,self.N)) + #self.Kmmi, Lm, Lmi, Kmm_logdet = pdinv(self.Kmm) + #self.dL_dpsi1 -= mdot(self.Kmmi,self.psi1*self.likelihood.precision.flatten().reshape(1,self.N)) #dB + + #########333333 + #self.Bi, self.LB, self.LBi, self.B_logdet = pdinv(self.B) + #########333333 + + + else: raise NotImplementedError, "homoscedastic fitc not implemented" # Remove extra term from dL_dpsi1 @@ -140,8 +158,11 @@ class generalized_FITC(sparse_GP): A = -0.5*self.N*self.D*np.log(2.*np.pi) +0.5*np.sum(np.log(self.likelihood.precision)) -0.5*np.sum(self.V*self.likelihood.Y) else: A = -0.5*self.N*self.D*(np.log(2.*np.pi) + np.log(self.likelihood._variance)) -0.5*self.likelihood.precision*self.likelihood.trYYT - C = -0.5*self.D * (self.B_logdet + self.M*np.log(sf2)) - D = 0.5*np.trace(self.Cpsi1VVpsi1) + C = -self.D * (np.sum(np.log(np.diag(self.LB))) + 0.5*self.M*np.log(sf2)) + #C = -0.5*self.D * (self.B_logdet + self.M*np.log(sf2)) + D = 0.5*np.sum(np.square(self._LBi_Lmi_psi1V)) + #self.Cpsi1VVpsi1 = np.dot(self.Cpsi1V,self.psi1V.T) + #D_ = 0.5*np.trace(self.Cpsi1VVpsi1) return A+C+D def _raw_predict(self, Xnew, which_parts, full_cov=False):