From b6bac53da815f7cd5e7920d3faac58cc1c9424a6 Mon Sep 17 00:00:00 2001 From: Neil Lawrence Date: Sun, 5 May 2013 07:59:07 +0100 Subject: [PATCH 1/6] Minor changes --- GPy/core/model.py | 7 ++++--- GPy/examples/dimensionality_reduction.py | 4 ++-- 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/GPy/core/model.py b/GPy/core/model.py index 25c10b42..ef05a2cb 100644 --- a/GPy/core/model.py +++ b/GPy/core/model.py @@ -416,9 +416,10 @@ class model(parameterised): """ return an array describing the sesitivity of the model to each input - NB. Right now, we're basing this on the lengthscales (or variances) of the kernel. - TODO: proper sensitivity analysis - """ + NB. Right now, we're basing this on the lengthscales (or + variances) of the kernel. TODO: proper sensitivity analysis + where we integrate across the model inputs and evaluate the + effect on the variance of the model output. """ if not hasattr(self, 'kern'): raise ValueError, "this model has no kernel" diff --git a/GPy/examples/dimensionality_reduction.py b/GPy/examples/dimensionality_reduction.py index 75820407..931e2eed 100644 --- a/GPy/examples/dimensionality_reduction.py +++ b/GPy/examples/dimensionality_reduction.py @@ -82,11 +82,11 @@ def BGPLVM_oil(optimize=True, N=100, Q=10, M=15, max_f_eval=300): m.ensure_default_constraints() y = m.likelihood.Y[0, :] - fig,(latent_axes,hist_axes) = plt.subplots(1,2) + fig,(latent_axes,sense_axes) = plt.subplots(1,2) plt.sca(latent_axes) m.plot_latent() data_show = GPy.util.visualize.vector_show(y) - lvm_visualizer = GPy.util.visualize.lvm_dimselect(m.X[0, :], m, data_show, latent_axes=latent_axes, hist_axes=hist_axes) + lvm_visualizer = GPy.util.visualize.lvm_dimselect(m.X[0, :], m, data_show, latent_axes=latent_axes, sense_axes=sense_axes) raw_input('Press enter to finish') plt.close('all') # # plot From 7f138b8b01973531f0e9b090861bc82828c8410e Mon Sep 17 00:00:00 2001 From: James Hensman Date: Tue, 7 May 2013 12:49:39 +0100 Subject: [PATCH 2/6] much tidy9ing in sparse_GP --- GPy/models/sparse_GP.py | 97 ++++++++++++++++------------------------- 1 file changed, 38 insertions(+), 59 deletions(-) diff --git a/GPy/models/sparse_GP.py b/GPy/models/sparse_GP.py index f04c9bd5..70f3899f 100644 --- a/GPy/models/sparse_GP.py +++ b/GPy/models/sparse_GP.py @@ -3,7 +3,7 @@ import numpy as np import pylab as pb -from ..util.linalg import mdot, jitchol, chol_inv, pdinv, trace_dot, tdot +from ..util.linalg import mdot, jitchol, tdot, symmetrify from ..util.plot import gpplot from .. import kern from GP import GP @@ -68,13 +68,11 @@ class sparse_GP(GP): self.psi2 = None def _computations(self): - #TODO: find routine to multiply triangular matrices - sf = self.scale_factor sf2 = sf**2 - #invert Kmm - self.Kmmi, self.Lm, self.Lmi, self.Kmm_logdet = pdinv(self.Kmm) + #factor Kmm + self.Lm = jitchol(self.Kmm) #The rather complex computations of self.A if self.likelihood.is_heteroscedastic: @@ -90,7 +88,6 @@ class sparse_GP(GP): self.A = tdot(tmp) else: tmp = self.psi1*(np.sqrt(self.likelihood.precision.flatten().reshape(1,self.N))/sf) - #self.psi2_beta_scaled = tdot(tmp) tmp, _ = linalg.lapack.flapack.dtrtrs(self.Lm,np.asfortranarray(tmp),lower=1) self.A = tdot(tmp) else: @@ -101,20 +98,16 @@ class sparse_GP(GP): if not np.allclose(evals, clipped_evals): print "Warning: clipping posterior eigenvalues" tmp = evecs*np.sqrt(clipped_evals) - #self.psi2_beta_scaled = tdot(tmp) tmp, _ = linalg.lapack.flapack.dtrtrs(self.Lm,np.asfortranarray(tmp),lower=1) self.A = tdot(tmp) else: tmp = self.psi1*(np.sqrt(self.likelihood.precision)/sf) - #self.psi2_beta_scaled = tdot(tmp) tmp, _ = linalg.lapack.flapack.dtrtrs(self.Lm,np.asfortranarray(tmp),lower=1) self.A = tdot(tmp) - #invert B and compute C. C is the posterior covariance of u + #factor B self.B = np.eye(self.M)/sf2 + self.A - self.Bi, self.LB, self.LBi, self.B_logdet = pdinv(self.B) - tmp = linalg.lapack.flapack.dtrtrs(self.Lm,np.asfortranarray(self.Bi),lower=1,trans=1)[0] - self.C = linalg.lapack.flapack.dtrtrs(self.Lm,np.asfortranarray(tmp.T),lower=1,trans=1)[0] + self.LB = jitchol(self.B) self.V = (self.likelihood.precision/self.scale_factor)*self.likelihood.Y self.psi1V = np.dot(self.psi1, self.V) @@ -122,41 +115,8 @@ class sparse_GP(GP): #back substutue C into psi1V tmp,info1 = linalg.lapack.flapack.dtrtrs(self.Lm,np.asfortranarray(self.psi1V),lower=1,trans=0) self._LBi_Lmi_psi1V,_ = linalg.lapack.flapack.dtrtrs(self.LB,np.asfortranarray(tmp),lower=1,trans=0) - self._P = tdot(tmp) tmp,info2 = linalg.lapack.flapack.dpotrs(self.LB,tmp,lower=1) self.Cpsi1V,info3 = linalg.lapack.flapack.dtrtrs(self.Lm,tmp,lower=1,trans=1) - #self.Cpsi1V = np.dot(self.C,self.psi1V) - - self.E = tdot(self.Cpsi1V/sf) - - - - # Compute dL_dpsi # FIXME: this is untested for the heterscedastic + uncertin inputs case - self.dL_dpsi0 = - 0.5 * self.D * (self.likelihood.precision * np.ones([self.N,1])).flatten() - self.dL_dpsi1 = np.dot(self.Cpsi1V,self.V.T) - if self.likelihood.is_heteroscedastic: - if self.has_uncertain_inputs: - #self.dL_dpsi2 = 0.5 * self.likelihood.precision[:,None,None] * self.D * self.Kmmi[None,:,:] # dB - #self.dL_dpsi2 += - 0.5 * self.likelihood.precision[:,None,None]/sf2 * self.D * self.C[None,:,:] # dC - #self.dL_dpsi2 += - 0.5 * self.likelihood.precision[:,None,None]* self.E[None,:,:] # dD - self.dL_dpsi2 = 0.5*self.likelihood.precision[:,None,None]*(self.D*(self.Kmmi - self.C/sf2) -self.E)[None,:,:] - else: - #self.dL_dpsi1 += mdot(self.Kmmi,self.psi1*self.likelihood.precision.flatten().reshape(1,self.N)) #dB - #self.dL_dpsi1 += -mdot(self.C,self.psi1*self.likelihood.precision.flatten().reshape(1,self.N)/sf2) #dC - #self.dL_dpsi1 += -mdot(self.E,self.psi1*self.likelihood.precision.flatten().reshape(1,self.N)) #dD - self.dL_dpsi1 += np.dot(self.Kmmi - self.C/sf2 -self.E,self.psi1*self.likelihood.precision.reshape(1,self.N)) - self.dL_dpsi2 = None - - else: - self.dL_dpsi2 = 0.5*self.likelihood.precision*(self.D*(self.Kmmi - self.C/sf2) -self.E) - if self.has_uncertain_inputs: - #repeat for each of the N psi_2 matrices - self.dL_dpsi2 = np.repeat(self.dL_dpsi2[None,:,:],self.N,axis=0) - else: - #subsume back into psi1 (==Kmn) - self.dL_dpsi1 += 2.*np.dot(self.dL_dpsi2,self.psi1) - self.dL_dpsi2 = None - # Compute dL_dKmm tmp = tdot(self._LBi_Lmi_psi1V) @@ -166,23 +126,38 @@ class sparse_GP(GP): tmp += self.D*np.eye(self.M) self.dL_dKmm = backsub_both_sides(self.Lm,tmp) + # Compute dL_dpsi # FIXME: this is untested for the heterscedastic + uncertain inputs case + self.dL_dpsi0 = - 0.5 * self.D * (self.likelihood.precision * np.ones([self.N,1])).flatten() + self.dL_dpsi1 = np.dot(self.Cpsi1V,self.V.T) + dL_dpsi2_beta = 0.5*backsub_both_sides(self.Lm,self.D*np.eye(self.M) - self.DBi_plus_BiPBi) + if self.likelihood.is_heteroscedastic: + if self.has_uncertain_inputs: + self.dL_dpsi2 = self.likelihood.precision[:,None,None]*dL_dpsi2_beta[None,:,:] + else: + self.dL_dpsi1 += 2.*np.dot(dL_dpsi2_beta,self.psi1*self.likelihood.precision.reshape(1,self.N)) + self.dL_dpsi2 = None + else: + dL_dpsi2 = self.likelihood.precision*dL_dpsi2_beta + if self.has_uncertain_inputs: + #repeat for each of the N psi_2 matrices + self.dL_dpsi2 = np.repeat(dL_dpsi2[None,:,:],self.N,axis=0) + else: + #subsume back into psi1 (==Kmn) + self.dL_dpsi1 += 2.*np.dot(dL_dpsi2,self.psi1) + self.dL_dpsi2 = None + + #the partial derivative vector for the likelihood if self.likelihood.Nparams ==0: #save computation here. self.partial_for_likelihood = None elif self.likelihood.is_heteroscedastic: raise NotImplementedError, "heteroscedatic derivates not implemented" - #self.partial_for_likelihood = - 0.5 * self.D*self.likelihood.precision + 0.5 * (self.likelihood.Y**2).sum(1)*self.likelihood.precision**2 #dA - #self.partial_for_likelihood += 0.5 * self.D * (self.psi0*self.likelihood.precision**2 - (self.psi2*self.Kmmi[None,:,:]*self.likelihood.precision[:,None,None]**2).sum(1).sum(1)/sf2) #dB - #self.partial_for_likelihood += 0.5 * self.D * np.sum(self.Bi*self.A)*self.likelihood.precision #dC - #self.partial_for_likelihood += -np.diag(np.dot((self.C - 0.5 * mdot(self.C,self.psi2_beta_scaled,self.C) ) , self.psi1VVpsi1 ))*self.likelihood.precision #dD else: #likelihood is not heterscedatic self.partial_for_likelihood = - 0.5 * self.N*self.D*self.likelihood.precision + 0.5 * self.likelihood.trYYT*self.likelihood.precision**2 self.partial_for_likelihood += 0.5 * self.D * (self.psi0.sum()*self.likelihood.precision**2 - np.trace(self.A)*self.likelihood.precision*sf2) - #self.partial_for_likelihood += 0.5 * self.D * trace_dot(self.Bi,self.A)*self.likelihood.precision - #self.partial_for_likelihood += self.likelihood.precision*(0.5*trace_dot(self.psi2_beta_scaled,self.E*sf2) - np.sum(np.square(self._LBi_Lmi_psi1V))) - self.partial_for_likelihood += self.likelihood.precision*(0.5*trace_dot(self.A,self.DBi_plus_BiPBi) - np.sum(np.square(self._LBi_Lmi_psi1V))) + self.partial_for_likelihood += self.likelihood.precision*(0.5*np.sum(self.A*self.DBi_plus_BiPBi) - np.sum(np.square(self._LBi_Lmi_psi1V))) @@ -195,7 +170,7 @@ class sparse_GP(GP): 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 B = -0.5*self.D*(np.sum(self.likelihood.precision*self.psi0) - np.trace(self.A)*sf2) - C = -0.5*self.D * (self.B_logdet + self.M*np.log(sf2)) + C = -self.D * (np.sum(np.log(np.diag(self.LB))) + 0.5*self.M*np.log(sf2)) D = 0.5*np.sum(np.square(self._LBi_Lmi_psi1V)) return A+B+C+D @@ -259,22 +234,26 @@ class sparse_GP(GP): """ dL_dZ = 2.*self.kern.dK_dX(self.dL_dKmm, self.Z) # factor of two becase of vertical and horizontal 'stripes' in dKmm_dZ if self.has_uncertain_inputs: - dL_dZ += self.kern.dpsi1_dZ(self.dL_dpsi1,self.Z,self.X, self.X_variance) + dL_dZ += self.kern.dpsi1_dZ(self.dL_dpsi1, self.Z, self.X, self.X_variance) dL_dZ += self.kern.dpsi2_dZ(self.dL_dpsi2, self.Z, self.X, self.X_variance) else: - dL_dZ += self.kern.dK_dX(self.dL_dpsi1,self.Z,self.X) + dL_dZ += self.kern.dK_dX(self.dL_dpsi1, self.Z, self.X) return dL_dZ def _raw_predict(self, Xnew, which_parts='all', full_cov=False): """Internal helper function for making predictions, does not account for normalization""" - Kx = self.kern.K(self.Z, Xnew) - mu = mdot(Kx.T, self.C/self.scale_factor, self.psi1V) + Bi,_ = linalg.lapack.flapack.dpotri(self.LB,lower=0) # WTH? this lower switch should be 1, but that doesn't work! + symmetrify(Bi) + Kmmi_LmiBLmi = backsub_both_sides(self.Lm,np.eye(self.M) - Bi) + + Kx = self.kern.K(self.Z, Xnew, which_parts=which_parts) + mu = np.dot(Kx.T, self.Cpsi1V/self.scale_factor) if full_cov: Kxx = self.kern.K(Xnew,which_parts=which_parts) - var = Kxx - mdot(Kx.T, (self.Kmmi - self.C/self.scale_factor**2), Kx) #NOTE this won't work for plotting + var = Kxx - mdot(Kx.T, Kmmi_LmiBLmi, Kx) #NOTE this won't work for plotting else: Kxx = self.kern.Kdiag(Xnew,which_parts=which_parts) - var = Kxx - np.sum(Kx*np.dot(self.Kmmi - self.C/self.scale_factor**2, Kx),0) + var = Kxx - np.sum(Kx*np.dot(Kmmi_LmiBLmi, Kx),0) return mu,var[:,None] From ce2884f0a7dc26087a5225bc92e39643920e3e16 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Tue, 7 May 2013 18:02:10 +0100 Subject: [PATCH 3/6] weaved linear kern --- GPy/kern/linear.py | 107 +++++++++++++++++++++++++++++++++++---------- 1 file changed, 83 insertions(+), 24 deletions(-) diff --git a/GPy/kern/linear.py b/GPy/kern/linear.py index 396b1aec..16ef2499 100644 --- a/GPy/kern/linear.py +++ b/GPy/kern/linear.py @@ -5,6 +5,7 @@ from kernpart import kernpart import numpy as np from ..util.linalg import tdot +from scipy import weave class linear(kernpart): """ @@ -171,33 +172,91 @@ class linear(kernpart): self._psi_computations(Z, mu, S) AZZA = self.ZA.T[:, None, :, None] * self.ZA[None, :, None, :] AZZA = AZZA + AZZA.swapaxes(1, 2) - target_S += (dL_dpsi2[:, :, :, None] * self.ZA[None, :, None, :] * self.ZA[None, None, :, :]).sum(1).sum(1) - dpsi2_dmu = (dL_dpsi2[:, :, :, None] * np.tensordot(mu, AZZA, (-1, 0))).sum(1).sum(1) - target_mu += dpsi2_dmu + AZZA_2 = AZZA/2. + #muAZZA = np.tensordot(mu,AZZA,(-1,0)) + #target_mu_dummy, target_S_dummy = np.zeros_like(target_mu), np.zeros_like(target_S) + #target_mu_dummy += (dL_dpsi2[:, :, :, None] * muAZZA).sum(1).sum(1) + #target_S_dummy += (dL_dpsi2[:, :, :, None] * self.ZA[None, :, None, :] * self.ZA[None, None, :, :]).sum(1).sum(1) + + #Using weave, we can exploiut the symmetry of this problem: + code = """ + int n, m, mm,q,qq; + double factor,tmp; + #pragma omp parallel for private(m,mm,q,qq,factor,tmp) + for(n=0;n + #include + """ + weave_options = {'headers' : [''], + 'extra_compile_args': ['-fopenmp -O3'], #-march=native'], + 'extra_link_args' : ['-lgomp']} + + N,M,Q = mu.shape[0],Z.shape[0],mu.shape[1] + weave.inline(code, support_code=support_code, libraries=['gomp'], + arg_names=['N','M','Q','mu','AZZA','AZZA_2','target_mu','target_S','dL_dpsi2'], + type_converters=weave.converters.blitz,**weave_options) + def dpsi2_dZ(self, dL_dpsi2, Z, mu, S, target): self._psi_computations(Z, mu, S) -# mu2_S = np.sum(self.mu2_S, 0) # Q, -# import ipdb;ipdb.set_trace() -# psi2_dZ_real = np.zeros((mu.shape[0], Z.shape[0], Z.shape[1])) -# for n in range(mu.shape[0]): -# for m in range(Z.shape[0]): -# tmp = self.variances * (tdot(self._mu[n:n + 1].T) + np.diag(S[n])) -# psi2_dZ_real[n, m, :] = np.dot(tmp, ( -# self._Z[m:m + 1] * self.variances).T).T -# tmp = self._Z[m:m + 1] * self.variances -# tmp = np.dot(tmp, (tdot(self._mu[n:n + 1].T) + np.diag(S[n]))) -# psi2_dZ_real[n, m, :] = tmp * self.variances -# for m_prime in range(Z.shape[0]): -# if m == m_prime: -# psi2_dZ_real[n, m, :] *= 2 -# prod = (dL_dpsi2[:, :, :, None] * np.eye(Z.shape[0])[None, :, :, None] * (self.ZAinner * self.variances).swapaxes(0, 1)[:, :, None, :]) -# psi2_dZ = prod.swapaxes(1, 2) + prod - psi2_dZ = dL_dpsi2[:, :, :, None] * self.variances * self.ZAinner[:, :, None, :] - target += psi2_dZ.sum(0).sum(0) -# import ipdb;ipdb.set_trace() -# psi2_dZ_old = (dL_dpsi2[:, :, :, None] * (self.mu2_S[:, None, None, :] * (Z * np.square(self.variances)[None, :])[None, None, :, :])).sum(0).sum(1) -# target += (dL_dpsi2[:, :, :, None] * psi2_dZ_real[:, :, None, :]).sum(0).sum(0) * 2 # (self.variances * np.dot(self.inner, self.ZA.T)).sum(1) + #psi2_dZ = dL_dpsi2[:, :, :, None] * self.variances * self.ZAinner[:, :, None, :] + #dummy_target = np.zeros_like(target) + #dummy_target += psi2_dZ.sum(0).sum(0) + + AZA = self.variances*self.ZAinner + code=""" + int n,m,mm,q; + #pragma omp parallel for private(n,mm,q) + for(m=0;m + #include + """ + weave_options = {'headers' : [''], + 'extra_compile_args': ['-fopenmp -O3'], #-march=native'], + 'extra_link_args' : ['-lgomp']} + + N,M,Q = mu.shape[0],Z.shape[0],mu.shape[1] + weave.inline(code, support_code=support_code, libraries=['gomp'], + arg_names=['N','M','Q','AZA','target','dL_dpsi2'], + type_converters=weave.converters.blitz,**weave_options) + + + + #---------------------------------------# # Precomputations # From 8f75f6d66b009f7c0a24274880c1df5adc3965df Mon Sep 17 00:00:00 2001 From: James Hensman Date: Tue, 7 May 2013 18:50:13 +0100 Subject: [PATCH 4/6] tried to eliminate the regexp overflow error for large models --- GPy/core/model.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/GPy/core/model.py b/GPy/core/model.py index 069c37b0..94202396 100644 --- a/GPy/core/model.py +++ b/GPy/core/model.py @@ -203,7 +203,7 @@ class model(parameterised): else: self._set_params_transformed(initial_parameters) - def ensure_default_constraints(self, warn=False): + def ensure_default_constraints(self): """ Ensure that any variables which should clearly be positive have been constrained somehow. """ @@ -214,11 +214,11 @@ class model(parameterised): for s in positive_strings: for i in self.grep_param_names(s): if not (i in currently_constrained): - to_make_positive.append(re.escape(param_names[i])) - if warn: - print "Warning! constraining %s positive" % s + #to_make_positive.append(re.escape(param_names[i])) + to_make_positive.append(i) if len(to_make_positive): - self.constrain_positive('(' + '|'.join(to_make_positive) + ')') + #self.constrain_positive('(' + '|'.join(to_make_positive) + ')') + self.constrain_positive(np.asarray(to_make_positive)) From b504c5f55084f478e85b182b6c826b14ec05eddf Mon Sep 17 00:00:00 2001 From: Ricardo Date: Tue, 7 May 2013 21:30:08 +0100 Subject: [PATCH 5/6] Broken file fixed --- GPy/likelihoods/EP.py | 12 +++--------- 1 file changed, 3 insertions(+), 9 deletions(-) diff --git a/GPy/likelihoods/EP.py b/GPy/likelihoods/EP.py index 8307b6b4..685195ba 100644 --- a/GPy/likelihoods/EP.py +++ b/GPy/likelihoods/EP.py @@ -196,9 +196,8 @@ class EP(likelihood): self.tau_tilde[i] = self.tau_tilde[i] + Delta_tau self.v_tilde[i] = self.v_tilde[i] + Delta_v #Posterior distribution parameters update - #LLT = LLT + np.outer(Kmn[:,i],Kmn[:,i])*Delta_tau - #L = jitchol(LLT) - cholupdate(L,Kmn[:,i]*np.sqrt(Delta_tau)) + LLT = LLT + np.outer(Kmn[:,i],Kmn[:,i])*Delta_tau + L = jitchol(LLT) V,info = linalg.lapack.flapack.dtrtrs(L,Kmn,lower=1) Sigma_diag = np.sum(V*V,-2) si = np.sum(V.T*V[:,i],-1) @@ -251,6 +250,7 @@ class EP(likelihood): R = R0.copy() Diag = Diag0.copy() Sigma_diag = Knn_diag + RPT0 = np.dot(R0,P0.T) """ Initial values - Cavity distribution parameters: @@ -306,13 +306,7 @@ class EP(likelihood): Iplus_Dprod_i = 1./(1.+ Diag0 * self.tau_tilde) Diag = Diag0 * Iplus_Dprod_i P = Iplus_Dprod_i[:,None] * P0 - - #Diag = Diag0/(1.+ Diag0 * self.tau_tilde) - #P = (Diag / Diag0)[:,None] * P0 - RPT0 = np.dot(R0,P0.T) L = jitchol(np.eye(M) + np.dot(RPT0,((1. - Iplus_Dprod_i)/Diag0)[:,None]*RPT0.T)) - #L = jitchol(np.eye(M) + np.dot(RPT0,(1./Diag0 - Iplus_Dprod_i/Diag0)[:,None]*RPT0.T)) - #L = jitchol(np.eye(M) + np.dot(RPT0,(1./Diag0 - Diag/(Diag0**2))[:,None]*RPT0.T)) R,info = linalg.lapack.flapack.dtrtrs(L,R0,lower=1) RPT = np.dot(R,P.T) Sigma_diag = Diag + np.sum(RPT.T*RPT.T,-1) From 71b845eb603833eba01ea80d5eaa3a0493011c9a Mon Sep 17 00:00:00 2001 From: Ricardo Date: Wed, 8 May 2013 07:09:00 +0100 Subject: [PATCH 6/6] 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):