apunkt/GPy
1
0
Fork 0
mirror of https://github.com/SheffieldML/GPy.git synced 2026-05-18 13:55:14 +02:00
Code Issues Projects Releases Packages Wiki Activity Actions

FITC is back

Browse source
This commit is contained in:
Ricardo 2013-06-05 16:30:41 +01:00
parent d238f7709f
commit 7fb3466fba
1 changed files with 32 additions and 98 deletions
Download patch file Download diff file Expand all files Collapse all files

130
GPy/core/fitc.py
View file

@ -7,14 +7,8 @@ from ..util.linalg import mdot, jitchol, chol_inv, tdot, symmetrify,pdinv
from ..util.plot import gpplot from ..util.plot import gpplot
from .. import kern from .. import kern
from scipy import stats, linalg from scipy import stats, linalg
from gp_base import GPBase
from sparse_gp import SparseGP 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): class FITC(SparseGP):
""" """
sparse FITC approximation sparse FITC approximation
@ -27,48 +21,13 @@ class FITC(SparseGP):
:type kernel: a GPy.kern.kern instance :type kernel: a GPy.kern.kern instance
:param Z: inducing inputs (optional, see note) :param Z: inducing inputs (optional, see note)
:type Z: np.ndarray (M x Q) | None :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) :param normalize_(X|Y) : whether to normalize the data before computing (predictions will be in original scales)
:type normalize_(X|Y): bool :type normalize_(X|Y): bool
""" """
def __init__(self, X, likelihood, kernel, Z, normalize_X=False): def __init__(self, X, likelihood, kernel, Z, normalize_X=False):
GPBase.__init__(self, X, likelihood, kernel, normalize_X=normalize_X) 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"
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()
def update_likelihood_approximation(self): def update_likelihood_approximation(self):
""" """
@ -77,47 +36,33 @@ class FITC(SparseGP):
For a Gaussian likelihood, no iteration is required: For a Gaussian likelihood, no iteration is required:
this function does nothing this function does nothing
""" """
if self.has_uncertain_inputs: self.likelihood.fit_FITC(self.Kmm,self.psi1,self.psi0)
raise NotImplementedError, "FITC approximation not implemented for uncertain inputs" self._set_params(self._get_params()) # update the GP
else:
self.likelihood.fit_FITC(self.Kmm,self.psi1,self.psi0)
self._set_params(self._get_params()) # update the GP
def _compute_kernel_matrices(self): def _compute_kernel_matrices(self):
# kernel computations, using BGPLVM notation # kernel computations, using BGPLVM notation
self.Kmm = self.kern.K(self.Z) self.Kmm = self.kern.K(self.Z)
if self.has_uncertain_inputs: self.psi0 = self.kern.Kdiag(self.X)
self.psi0 = self.kern.psi0(self.Z, self.X, self.X_variance) self.psi1 = self.kern.K(self.Z, self.X)
self.psi1 = self.kern.psi1(self.Z, self.X, self.X_variance).T self.psi2 = None
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
def _computations(self): def _computations(self):
#factor Kmm #factor Kmm
self.Lm = jitchol(self.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) Lmipsi1 = np.dot(self.Lmi,self.psi1)
self.Qnn = np.dot(Lmipsi1.T,Lmipsi1).copy() self.Qnn = np.dot(Lmipsi1.T,Lmipsi1).copy()
self.Diag0 = self.psi0 - np.diag(self.Qnn) 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 self.V_star = self.beta_star * self.likelihood.Y
# The rather complex computations of self.A # The rather complex computations of self.A
if self.has_uncertain_inputs: tmp = self.psi1 * (np.sqrt(self.beta_star.flatten().reshape(1, self.N)))
raise NotImplementedError tmp, _ = linalg.lapack.flapack.dtrtrs(self.Lm, np.asfortranarray(tmp), lower=1)
else: self.A = tdot(tmp)
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)
# factor B # 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.LB = jitchol(self.B)
self.LBi = chol_inv(self.LB) self.LBi = chol_inv(self.LB)
self.psi1V = np.dot(self.psi1, self.V_star) 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): 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) 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),:] _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 _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._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_dtheta += self.kern.dK_dtheta(_dKmm,self.Z)
self._dKmm_dX += 2.*self.kern.dK_dX(_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,:]) 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. # save computation here.
self.partial_for_likelihood = None self.partial_for_likelihood = None
elif self.likelihood.is_heteroscedastic: elif self.likelihood.is_heteroscedastic:
raise NotImplementedError, "heteroscedatic derivates not implemented" raise NotImplementedError, "heteroscedatic derivates not implemented."
else: else:
# likelihood is not heterscedatic # likelihood is not heterscedatic
dbstar_dnoise = self.likelihood.precision * (self.beta_star**2 * self.Diag0[:,None] - self.beta_star) 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))) return np.hstack((self.dL_dZ().flatten(), self.dL_dtheta(), self.likelihood._gradients(partial=self.partial_for_likelihood)))
def dL_dtheta(self): def dL_dtheta(self):
if self.has_uncertain_inputs: dL_dtheta = self.kern.dKdiag_dtheta(self._dL_dpsi0,self.X)
raise NotImplementedError, "FITC approximation not implemented for uncertain inputs" dL_dtheta += self.kern.dK_dtheta(self._dL_dpsi1,self.X,self.Z)
else: dL_dtheta += self.kern.dK_dtheta(self._dL_dKmm,X=self.Z)
dL_dtheta = self.kern.dKdiag_dtheta(self._dL_dpsi0,self.X) dL_dtheta += self._dKmm_dtheta
dL_dtheta += self.kern.dK_dtheta(self._dL_dpsi1,self.X,self.Z) dL_dtheta += self._dpsi1_dtheta
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: dL_dZ = self.kern.dK_dX(self._dL_dpsi1.T,self.Z,self.X)
raise NotImplementedError, "FITC approximation not implemented for uncertain inputs" dL_dZ += 2. * self.kern.dK_dX(self._dL_dKmm,X=self.Z)
else: dL_dZ += self._dpsi1_dX
dL_dZ = self.kern.dK_dX(self._dL_dpsi1.T,self.Z,self.X) dL_dZ += self._dKmm_dX
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, 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: if self.likelihood.is_heteroscedastic:
Iplus_Dprod_i = 1./(1.+ self.Diag0 * self.likelihood.precision.flatten()) Iplus_Dprod_i = 1./(1.+ self.Diag0 * self.likelihood.precision.flatten())
self.Diag = self.Diag0 * Iplus_Dprod_i self.Diag = self.Diag0 * Iplus_Dprod_i
self.P = Iplus_Dprod_i[:,None] * self.psi1.T self.P = Iplus_Dprod_i[:,None] * self.psi1.T
self.RPT0 = np.dot(self.Lmi,self.psi1) 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.num_inducing) + 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.R,info = linalg.lapack.flapack.dtrtrs(self.L,self.Lmi,lower=1)
self.RPT = np.dot(self.R,self.P.T) self.RPT = np.dot(self.R,self.P.T)
self.Sigma = np.diag(self.Diag) + np.dot(self.RPT.T,self.RPT) self.Sigma = np.diag(self.Diag) + np.dot(self.RPT.T,self.RPT)
self.w = self.Diag * self.likelihood.v_tilde self.w = self.Diag * self.likelihood.v_tilde
@ -275,8 +209,8 @@ class FITC(SparseGP):
# = I - [RPT0] * (U*U.T)^-1 * [RPT0].T # = I - [RPT0] * (U*U.T)^-1 * [RPT0].T
# = I - V.T * V # = I - V.T * V
U = np.linalg.cholesky(np.diag(self.Diag0) + self.Qnn) U = np.linalg.cholesky(np.diag(self.Diag0) + self.Qnn)
V,info = linalg.flapack.dtrtrs(U,self.RPT0.T,lower=1) V,info = linalg.lapack.flapack.dtrtrs(U,self.RPT0.T,lower=1)
C = np.eye(self.M) - np.dot(V.T,V) C = np.eye(self.num_inducing) - np.dot(V.T,V)
mu_u = np.dot(C,self.RPT0)*(1./self.Diag0[None,:]) mu_u = np.dot(C,self.RPT0)*(1./self.Diag0[None,:])
#self.C = C #self.C = C
#self.RPT0 = np.dot(self.R0,self.Knm.T) P0.T #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) mu_star = np.dot(KR0T,mu_H)
if full_cov: if full_cov:
Kxx = self.kern.K(Xnew,which_parts=which_parts) 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: else:
Kxx = self.kern.Kdiag(Xnew,which_parts=which_parts) 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 return mu_star[:,None],var
else: else:
raise NotImplementedError, "homoscedastic FITC not implemented" raise NotImplementedError, "Heteroscedastic case not implemented."
""" """
Kx = self.kern.K(self.Z, Xnew) Kx = self.kern.K(self.Z, Xnew)
mu = mdot(Kx.T, self.C/self.scale_factor, self.psi1V) mu = mdot(Kx.T, self.C/self.scale_factor, self.psi1V)