Some changes

This commit is contained in:
Ricardo 2013-05-15 02:01:08 +01:00
parent 3cbe492678
commit 4e827e913d

View file

@ -62,8 +62,9 @@ class FITC(sparse_GP):
self.psi1V = np.dot(self.psi1, self.V_star) self.psi1V = np.dot(self.psi1, self.V_star)
# back substutue C into psi1V # back substutue C into psi1V
tmp, info1 = linalg.lapack.flapack.dtrtrs(self.Lm, np.asfortranarray(self.psi1V), lower=1, trans=0) Lmi_psi1V, 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._LBi_Lmi_psi1V, _ = linalg.lapack.flapack.dtrtrs(self.LB, np.asfortranarray(Lmi_psi1V), lower=1, trans=0)
Kmmipsi1 = np.dot(self.Lmi.T,Lmipsi1) Kmmipsi1 = np.dot(self.Lmi.T,Lmipsi1)
b_psi1_Ki = self.beta_star * Kmmipsi1.T b_psi1_Ki = self.beta_star * Kmmipsi1.T
@ -76,100 +77,50 @@ class FITC(sparse_GP):
VV_p_Ki = np.dot(VVT,Kmmipsi1.T) VV_p_Ki = np.dot(VVT,Kmmipsi1.T)
Ki_pVVp_Ki = np.dot(Kmmipsi1,VV_p_Ki) Ki_pVVp_Ki = np.dot(Kmmipsi1,VV_p_Ki)
psi1beta = self.psi1*self.beta_star.T psi1beta = self.psi1*self.beta_star.T
H = self.Kmm + mdot(self.psi1,self.beta_star*self.psi1.T) H = self.Kmm + mdot(self.psi1,psi1beta.T)
Hi, LH, LHi, logdetH = pdinv(H) Hi, LH, LHi, logdetH = pdinv(H)
betapsi1TLmiLBi = np.dot(psi1beta.T,LBiLmi.T) betapsi1TLmiLBi = np.dot(psi1beta.T,LBiLmi.T)
alpha = np.array([np.dot(a.T,a) for a in betapsi1TLmiLBi])[:,None] alpha = np.array([np.dot(a.T,a) for a in betapsi1TLmiLBi])[:,None]
gamma_1 = mdot(VVT,self.psi1.T,Hi) gamma_1 = mdot(VVT,self.psi1.T,Hi)
pHip = mdot(self.psi1.T,Hi,self.psi1) pHip = mdot(self.psi1.T,Hi,self.psi1)
gamma_2 = mdot(self.beta_star*pHip,self.V_star) gamma_2 = mdot(self.beta_star*pHip,self.V_star)
gamma_3 = self.V_star * mdot(self.V_star.T,pHip*self.beta_star).T #gamma_3 = self.V_star * mdot(self.V_star.T,pHip*self.beta_star).T
gamma_3 = self.V_star * gamma_2
self._dL_dpsi0 = -0.5 * self.beta_star#dA_dpsi0: logdet(self.beta_star)
self._dL_dpsi0 = .5 * self.V_star**2 #dA_psi0???? self._dL_dpsi0 += .5 * self.V_star**2 #dA_psi0: yT*beta_star*y
self._dL_dpsi0 += -0.5 * self.beta_star#dA_dpsi0
self._dL_dpsi0 += .5 *alpha #dC_dpsi0 self._dL_dpsi0 += .5 *alpha #dC_dpsi0
self._dL_dpsi0 += 0.5*mdot(self.beta_star*pHip,self.V_star)**2 - self.V_star * mdot(self.V_star.T,pHip*self.beta_star).T #dD_dpsi0 self._dL_dpsi0 += 0.5*mdot(self.beta_star*pHip,self.V_star)**2 - self.V_star * mdot(self.V_star.T,pHip*self.beta_star).T #dD_dpsi0
self._dL_dpsi1 = b_psi1_Ki.copy() #dA_dpsi1 self._dL_dpsi1 = b_psi1_Ki.copy() #dA_dpsi1: logdet(self.beta_star)
self._dL_dpsi1 += -np.dot(psi1beta.T,LBL_inv) #dC_dpsi1 self._dL_dpsi1 += -np.dot(psi1beta.T,LBL_inv) #dC_dpsi1
self._dL_dpsi1 += gamma_1 - mdot(psi1beta.T,Hi,self.psi1,gamma_1) #dD_dpsi1 self._dL_dpsi1 += gamma_1 - mdot(psi1beta.T,Hi,self.psi1,gamma_1) #dD_dpsi1
self._dL_dKmm = -0.5 * np.dot(Kmmipsi1,b_psi1_Ki) #dA_dKmm self._dL_dKmm = -0.5 * np.dot(Kmmipsi1,b_psi1_Ki) #dA_dKmm: logdet(self.beta_star)
self._dL_dKmm += -.5*Kmmi + .5*LBL_inv + mdot(LBL_inv,psi1beta,Kmmipsi1.T) #dC_dKmm self._dL_dKmm += -.5*Kmmi + .5*LBL_inv + mdot(LBL_inv,psi1beta,Kmmipsi1.T) #dC_dKmm
self._dL_dKmm += -.5 * mdot(Hi,self.psi1,gamma_1) #dD_dKmm self._dL_dKmm += -.5 * mdot(Hi,self.psi1,gamma_1) #dD_dKmm
_dpsi1_dtheta = 0 self._dpsi1_dtheta = 0
_dpsi1_dX = 0 self._dpsi1_dX = 0
_dKmm_dtheta = 0 self._dKmm_dtheta = 0
_dKmm_dX = 0 self._dKmm_dX = 0
for psi1_n,V_n,X_n,alpha_n,gamma_n,gamma_k in zip(self.psi1.T,self.V_star,self.X,alpha,gamma_2,gamma_3): for psi1_n,V_n,X_n,alpha_n,gamma_n,gamma_k in zip(self.psi1.T,self.V_star,self.X,alpha,gamma_2,gamma_3):
psin_K = np.dot(psi1_n[None,:],Kmmi) psin_K = np.dot(psi1_n[None,:],Kmmi)
_dpsi1 = -V_n**2 * psin_K #Diag_dA_dpsi1 _dpsi1 = -V_n**2 * psin_K #dA_dpsi1: yT*beta_star*y
_dpsi1 += - alpha_n * psin_K #Diag_dC_dpsi1 _dpsi1 += - alpha_n * psin_K #Diag_dC_dpsi1
_dpsi1 += - gamma_n**2 * psin_K + 2. * gamma_k * psin_K #Diag_dD_dpsi1 _dpsi1 += - gamma_n**2 * psin_K + 2. * gamma_k * psin_K #Diag_dD_dpsi1
_dKmm = .5*V_n**2 * np.dot(psin_K.T,psin_K) #Diag_dA_dKmm _dKmm = .5*V_n**2 * np.dot(psin_K.T,psin_K) #dA_dKmm: yT*beta_star*y
_dKmm += .5 * alpha_n * np.dot(psin_K.T,psin_K) #Diag_dC_dKmm _dKmm += .5 * alpha_n * np.dot(psin_K.T,psin_K) #Diag_dC_dKmm
_dKmm += .5*gamma_n**2 * np.dot(psin_K.T,psin_K) - gamma_n * np.dot(psin_K.T,psin_K) #Diag_dD_dKmm _dKmm += .5*gamma_n**2 * np.dot(psin_K.T,psin_K) - gamma_k * np.dot(psin_K.T,psin_K) #Diag_dD_dKmm
_dpsi1_dtheta += self.kern.dK_dtheta(_dpsi1,X_n[None,:],self.Z) self._dpsi1_dtheta += self.kern.dK_dtheta(_dpsi1,X_n[None,:],self.Z)
_dKmm_dtheta += self.kern.dK_dtheta(_dKmm,self.Z) self._dKmm_dtheta += self.kern.dK_dtheta(_dKmm,self.Z)
_dKmm_dX += 2.*self.kern.dK_dX(_dKmm ,self.Z)
_dpsi1_dX += self.kern.dK_dX(_dpsi1.T,self.Z,X_n[None,:])
#_dA_dpsi1 = -V_n**2 * np.dot(psi1_n[None,:],Kmmi)
#_dC_dpsi1 = - alpha_n * np.dot(psi1_n[None,:],Kmmi)
#_dD_dpsi1_1 = - gamma_n**2 * np.dot(psi1_n[None,:],Kmmi)
#_dD_dpsi1_2 = 2. * gamma_k * np.dot(psi1_n[None,:],Kmmi)
#dA_dpsi1_theta += self.kern.dK_dtheta(_dA_dpsi1,X_n[None,:],self.Z)
#_dC_dpsi1_dtheta += self.kern.dK_dtheta(_dC_dpsi1,X_n[None,:],self.Z)
#_dD_dpsi1_dtheta_1 += self.kern.dK_dtheta(_dD_dpsi1_1,X_n[None,:],self.Z)
#_dD_dpsi1_dtheta_2 += self.kern.dK_dtheta(_dD_dpsi1_2,X_n[None,:],self.Z)
#dA_dpsi1_X += self.kern.dK_dX(_dA_dpsi1.T,self.Z,X_n[None,:])
#_dC_dpsi1_dX += self.kern.dK_dX(_dC_dpsi1.T,self.Z,X_n[None,:])
#_dD_dpsi1_dX_1 += self.kern.dK_dX(_dD_dpsi1_1.T,self.Z,X_n[None,:])
#_dD_dpsi1_dX_2 += self.kern.dK_dX(_dD_dpsi1_2.T,self.Z,X_n[None,:])
#_dA_dKmm = .5*V_n**2 * np.dot(psin_K.T,psin_K)
#_dC_dKmm = .5 * alpha_n * np.dot(psin_K.T,psin_K)
#_dD_dKmm_1 = .5*gamma_n**2 * np.dot(psin_K.T,psin_K)
#_dD_dKmm_2 = - gamma_n * np.dot(psin_K.T,psin_K)
#dA_dKmm_theta += self.kern.dK_dtheta(_dA_dKmm,self.Z)
#_dC_dKmm_dtheta += self.kern.dK_dtheta(_dC_dKmm,self.Z)
#_dD_dKmm_dtheta_1 += self.kern.dK_dtheta(_dD_dKmm_1,self.Z)
#_dD_dKmm_dtheta_2 += self.kern.dK_dtheta(_dD_dKmm_2,self.Z)
#dA_dKmm_X += 2.*self.kern.dK_dX(_dA_dKmm,self.Z)
#_dC_dKmm_dX += 2.*self.kern.dK_dX(_dC_dKmm,self.Z)
#_dD_dKmm_dX_1 += 2.*self.kern.dK_dX(_dD_dKmm_1,self.Z)
#_dD_dKmm_dX_2 += 2.*self.kern.dK_dX(_dD_dKmm_2,self.Z)
self._dL_dtheta = self.kern.dKdiag_dtheta(self._dL_dpsi0,self.X)
self._dL_dtheta += self.kern.dK_dtheta(self._dL_dpsi1,self.X,self.Z)
self._dL_dtheta += self.kern.dK_dtheta(self._dL_dKmm,X=self.Z)
self._dL_dtheta += _dKmm_dtheta
self._dL_dtheta += _dpsi1_dtheta
self._dL_dX = self.kern.dK_dX(self._dL_dpsi1.T,self.Z,self.X)
self._dL_dX += 2. * self.kern.dK_dX(self._dL_dKmm,X=self.Z)
self._dL_dX += _dpsi1_dX
self._dL_dX += _dKmm_dX
self._dKmm_dX += 2.*self.kern.dK_dX(_dKmm ,self.Z)
self._dpsi1_dX += self.kern.dK_dX(_dpsi1.T,self.Z,X_n[None,:])
# the partial derivative vector for the likelihood # the partial derivative vector for the likelihood
@ -217,16 +168,21 @@ class FITC(sparse_GP):
if self.has_uncertain_inputs: if self.has_uncertain_inputs:
raise NotImplementedError, "FITC approximation not implemented for uncertain inputs" raise NotImplementedError, "FITC approximation not implemented for uncertain inputs"
else: else:
#dL_dtheta = self.dA_dtheta + self.dlogB_dtheta + self.dD_dtheta dL_dtheta = self.kern.dKdiag_dtheta(self._dL_dpsi0,self.X)
dL_dtheta = self._dL_dtheta #FIXME 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 return dL_dtheta
def dL_dZ(self): def dL_dZ(self):
if self.has_uncertain_inputs: if self.has_uncertain_inputs:
raise NotImplementedError, "FITC approximation not implemented for uncertain inputs" raise NotImplementedError, "FITC approximation not implemented for uncertain inputs"
else: else:
#dL_dZ = self.dA_dX + self.dlogB_dX + self.dD_dX dL_dZ = self.kern.dK_dX(self._dL_dpsi1.T,self.Z,self.X)
dL_dZ = self._dL_dX #FIXME 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 return dL_dZ
def _raw_predict(self, Xnew, which_parts, full_cov=False): def _raw_predict(self, Xnew, which_parts, full_cov=False):