From 7fb3466fbaf49bd105dbf500be98bdad75a33b76 Mon Sep 17 00:00:00 2001 From: Ricardo Date: Wed, 5 Jun 2013 16:30:41 +0100 Subject: [PATCH] FITC is back --- GPy/core/fitc.py | 130 ++++++++++++----------------------------------- 1 file changed, 32 insertions(+), 98 deletions(-) diff --git a/GPy/core/fitc.py b/GPy/core/fitc.py index a6bfb3d5..639e0561 100644 --- a/GPy/core/fitc.py +++ b/GPy/core/fitc.py @@ -7,14 +7,8 @@ from ..util.linalg import mdot, jitchol, chol_inv, tdot, symmetrify,pdinv from ..util.plot import gpplot from .. import kern from scipy import stats, linalg -from gp_base import GPBase from sparse_gp import SparseGP -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 FITC(SparseGP): """ sparse FITC approximation @@ -27,48 +21,13 @@ class FITC(SparseGP): :type kernel: a GPy.kern.kern instance :param Z: inducing inputs (optional, see note) :type Z: np.ndarray (M x Q) | None - :param M : Number of inducing points (optional, default 10. Ignored if Z is not None) - :type M: int :param normalize_(X|Y) : whether to normalize the data before computing (predictions will be in original scales) :type normalize_(X|Y): bool """ def __init__(self, X, likelihood, kernel, Z, normalize_X=False): - GPBase.__init__(self, X, likelihood, kernel, normalize_X=normalize_X) - - self.Z = Z - self.M = Z.shape[0] - self.likelihood = likelihood - - X_variance = None - if X_variance is None: - self.has_uncertain_inputs = False - else: - assert X_variance.shape == X.shape - self.has_uncertain_inputs = True - self.X_variance = X_variance - - if normalize_X: - self.Z = (self.Z.copy() - self._Xmean) / self._Xstd - - # normalize X uncertainty also - if self.has_uncertain_inputs: - self.X_variance /= np.square(self._Xstd) - - def _set_params(self, p): - self.Z = p[:self.M * self.input_dim].reshape(self.M, self.input_dim) - self.kern._set_params(p[self.Z.size:self.Z.size + self.kern.Nparam]) - self.likelihood._set_params(p[self.Z.size + self.kern.Nparam:]) - self._compute_kernel_matrices() - self._computations() - - def _get_params(self): - return np.hstack([self.Z.flatten(), self.kern._get_params_transformed(), self.likelihood._get_params()]) - - def _get_param_names(self): - return sum([['iip_%i_%i' % (i, j) for j in range(self.Z.shape[1])] for i in range(self.Z.shape[0])],[])\ - + self.kern._get_param_names_transformed() + self.likelihood._get_param_names() - + SparseGP.__init__(self, X, likelihood, kernel, Z, X_variance=None, normalize_X=False) + assert self.output_dim == 1, "FITC model is not defined for handling multiple outputs" def update_likelihood_approximation(self): """ @@ -77,47 +36,33 @@ class FITC(SparseGP): For a Gaussian likelihood, no iteration is required: this function does nothing """ - 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._set_params(self._get_params()) # update the GP + self.likelihood.fit_FITC(self.Kmm,self.psi1,self.psi0) + self._set_params(self._get_params()) # update the GP def _compute_kernel_matrices(self): # kernel computations, using BGPLVM notation self.Kmm = self.kern.K(self.Z) - if self.has_uncertain_inputs: - self.psi0 = self.kern.psi0(self.Z, self.X, self.X_variance) - self.psi1 = self.kern.psi1(self.Z, self.X, self.X_variance).T - self.psi2 = self.kern.psi2(self.Z, self.X, self.X_variance) - else: - self.psi0 = self.kern.Kdiag(self.X) - self.psi1 = self.kern.K(self.Z, self.X) - self.psi2 = None - + self.psi0 = self.kern.Kdiag(self.X) + self.psi1 = self.kern.K(self.Z, self.X) + self.psi2 = None def _computations(self): #factor Kmm self.Lm = jitchol(self.Kmm) - self.Lmi,info = linalg.lapack.flapack.dtrtrs(self.Lm,np.eye(self.M),lower=1) + self.Lmi,info = linalg.lapack.flapack.dtrtrs(self.Lm,np.eye(self.num_inducing),lower=1) Lmipsi1 = np.dot(self.Lmi,self.psi1) self.Qnn = np.dot(Lmipsi1.T,Lmipsi1).copy() self.Diag0 = self.psi0 - np.diag(self.Qnn) - self.beta_star = self.likelihood.precision/(1. + self.likelihood.precision*self.Diag0[:,None]) #Includes Diag0 in the precision + self.beta_star = self.likelihood.precision/(1. + self.likelihood.precision*self.Diag0[:,None]) #NOTE: beta_star contains Diag0 and the precision self.V_star = self.beta_star * self.likelihood.Y # The rather complex computations of self.A - if self.has_uncertain_inputs: - raise NotImplementedError - else: - if self.likelihood.is_heteroscedastic: - assert self.likelihood.output_dim == 1 - tmp = self.psi1 * (np.sqrt(self.beta_star.flatten().reshape(1, self.N))) - tmp, _ = linalg.lapack.flapack.dtrtrs(self.Lm, np.asfortranarray(tmp), lower=1) - self.A = tdot(tmp) + tmp = self.psi1 * (np.sqrt(self.beta_star.flatten().reshape(1, self.N))) + tmp, _ = linalg.lapack.flapack.dtrtrs(self.Lm, np.asfortranarray(tmp), lower=1) + self.A = tdot(tmp) # factor B - self.B = np.eye(self.M) + self.A + self.B = np.eye(self.num_inducing) + self.A self.LB = jitchol(self.B) self.LBi = chol_inv(self.LB) self.psi1V = np.dot(self.psi1, self.V_star) @@ -170,16 +115,10 @@ class FITC(SparseGP): for i,V_n,alpha_n,gamma_n,gamma_k in zip(range(self.N),self.V_star,alpha,gamma_2,gamma_3): K_pp_K = np.dot(Kmmipsi1[:,i:(i+1)],Kmmipsi1[:,i:(i+1)].T) - - #Diag_dpsi1 = Diag_dA_dpsi1: yT*beta_star*y + Diag_dC_dpsi1 +Diag_dD_dpsi1 _dpsi1 = (-V_n**2 - alpha_n + 2.*gamma_k - gamma_n**2) * Kmmipsi1.T[i:(i+1),:] - - #Diag_dKmm = Diag_dA_dKmm: yT*beta_star*y +Diag_dC_dKmm +Diag_dD_dKmm _dKmm = .5*(V_n**2 + alpha_n + gamma_n**2 - 2.*gamma_k) * K_pp_K #Diag_dD_dKmm - self._dpsi1_dtheta += self.kern.dK_dtheta(_dpsi1,self.X[i:i+1,:],self.Z) self._dKmm_dtheta += self.kern.dK_dtheta(_dKmm,self.Z) - self._dKmm_dX += 2.*self.kern.dK_dX(_dKmm ,self.Z) self._dpsi1_dX += self.kern.dK_dX(_dpsi1.T,self.Z,self.X[i:i+1,:]) @@ -188,7 +127,7 @@ class FITC(SparseGP): # save computation here. self.partial_for_likelihood = None elif self.likelihood.is_heteroscedastic: - raise NotImplementedError, "heteroscedatic derivates not implemented" + raise NotImplementedError, "heteroscedatic derivates not implemented." else: # likelihood is not heterscedatic dbstar_dnoise = self.likelihood.precision * (self.beta_star**2 * self.Diag0[:,None] - self.beta_star) @@ -225,35 +164,30 @@ class FITC(SparseGP): return np.hstack((self.dL_dZ().flatten(), self.dL_dtheta(), self.likelihood._gradients(partial=self.partial_for_likelihood))) def dL_dtheta(self): - if self.has_uncertain_inputs: - raise NotImplementedError, "FITC approximation not implemented for uncertain inputs" - else: - dL_dtheta = self.kern.dKdiag_dtheta(self._dL_dpsi0,self.X) - dL_dtheta += self.kern.dK_dtheta(self._dL_dpsi1,self.X,self.Z) - dL_dtheta += self.kern.dK_dtheta(self._dL_dKmm,X=self.Z) - dL_dtheta += self._dKmm_dtheta - dL_dtheta += self._dpsi1_dtheta + dL_dtheta = self.kern.dKdiag_dtheta(self._dL_dpsi0,self.X) + dL_dtheta += self.kern.dK_dtheta(self._dL_dpsi1,self.X,self.Z) + dL_dtheta += self.kern.dK_dtheta(self._dL_dKmm,X=self.Z) + dL_dtheta += self._dKmm_dtheta + dL_dtheta += self._dpsi1_dtheta return dL_dtheta def dL_dZ(self): - if self.has_uncertain_inputs: - raise NotImplementedError, "FITC approximation not implemented for uncertain inputs" - else: - dL_dZ = self.kern.dK_dX(self._dL_dpsi1.T,self.Z,self.X) - dL_dZ += 2. * self.kern.dK_dX(self._dL_dKmm,X=self.Z) - dL_dZ += self._dpsi1_dX - dL_dZ += self._dKmm_dX + dL_dZ = self.kern.dK_dX(self._dL_dpsi1.T,self.Z,self.X) + dL_dZ += 2. * self.kern.dK_dX(self._dL_dKmm,X=self.Z) + dL_dZ += self._dpsi1_dX + dL_dZ += self._dKmm_dX return dL_dZ - def _raw_predict(self, Xnew, which_parts, full_cov=False): + def _raw_predict(self, Xnew, X_variance_new=None, which_parts='all', full_cov=False): + assert X_variance_new is None, "FITC model is not defined for handling uncertain inputs." if self.likelihood.is_heteroscedastic: Iplus_Dprod_i = 1./(1.+ self.Diag0 * self.likelihood.precision.flatten()) self.Diag = self.Diag0 * Iplus_Dprod_i self.P = Iplus_Dprod_i[:,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.R,info = linalg.flapack.dtrtrs(self.L,self.Lmi,lower=1) + self.L = np.linalg.cholesky(np.eye(self.num_inducing) + np.dot(self.RPT0,((1. - Iplus_Dprod_i)/self.Diag0)[:,None]*self.RPT0.T)) + self.R,info = linalg.lapack.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) self.w = self.Diag * self.likelihood.v_tilde @@ -275,8 +209,8 @@ class FITC(SparseGP): # = I - [RPT0] * (U*U.T)^-1 * [RPT0].T # = I - V.T * V U = np.linalg.cholesky(np.diag(self.Diag0) + self.Qnn) - V,info = linalg.flapack.dtrtrs(U,self.RPT0.T,lower=1) - C = np.eye(self.M) - np.dot(V.T,V) + V,info = linalg.lapack.flapack.dtrtrs(U,self.RPT0.T,lower=1) + C = np.eye(self.num_inducing) - np.dot(V.T,V) mu_u = np.dot(C,self.RPT0)*(1./self.Diag0[None,:]) #self.C = C #self.RPT0 = np.dot(self.R0,self.Knm.T) P0.T @@ -292,13 +226,13 @@ class FITC(SparseGP): mu_star = np.dot(KR0T,mu_H) if full_cov: Kxx = self.kern.K(Xnew,which_parts=which_parts) - var = Kxx + np.dot(KR0T,np.dot(Sigma_H - np.eye(self.M),KR0T.T)) + var = Kxx + np.dot(KR0T,np.dot(Sigma_H - np.eye(self.num_inducing),KR0T.T)) else: Kxx = self.kern.Kdiag(Xnew,which_parts=which_parts) - var = (Kxx + np.sum(KR0T.T*np.dot(Sigma_H - np.eye(self.M),KR0T.T),0))[:,None] + var = (Kxx + np.sum(KR0T.T*np.dot(Sigma_H - np.eye(self.num_inducing),KR0T.T),0))[:,None] return mu_star[:,None],var else: - raise NotImplementedError, "homoscedastic FITC not implemented" + raise NotImplementedError, "Heteroscedastic case not implemented." """ Kx = self.kern.K(self.Z, Xnew) mu = mdot(Kx.T, self.C/self.scale_factor, self.psi1V)