more tidying in sparse_GP

This commit is contained in:
James Hensman 2013-05-20 09:38:08 +01:00
parent 3f22e61d2d
commit 025f31a43f

View file

@ -16,9 +16,9 @@ class sparse_GP(GP):
:param X: inputs
:type X: np.ndarray (N x Q)
:param likelihood: a likelihood instance, containing the observed data
:type likelihood: GPy.likelihood.(Gaussian | EP)
:param kernel : the kernel/covariance function. See link kernels
:type kernel: a GPy kernel
:type likelihood: GPy.likelihood.(Gaussian | EP | Laplace)
:param kernel : the kernel (covariance function). See link kernels
:type kernel: a GPy.kern.kern instance
:param X_variance: The uncertainty in the measurements of X (Gaussian variance)
:type X_variance: np.ndarray (N x Q) | None
:param Z: inducing inputs (optional, see note)
@ -30,8 +30,6 @@ class sparse_GP(GP):
"""
def __init__(self, X, likelihood, kernel, Z, X_variance=None, normalize_X=False):
# self.scale_factor = 100.0 # a scaling factor to help keep the algorithm stable
# self.auto_scale_factor = False
self.Z = Z
self.M = Z.shape[0]
self.likelihood = likelihood
@ -63,8 +61,6 @@ class sparse_GP(GP):
self.psi2 = None
def _computations(self):
# sf = self.scale_factor
# sf2 = sf ** 2
# factor Kmm
self.Lm = jitchol(self.Kmm)
@ -88,7 +84,6 @@ class sparse_GP(GP):
# factor B
# self.B = np.eye(self.M) / sf2 + self.A
self.B = np.eye(self.M) + self.A
self.LB = jitchol(self.B)
@ -104,8 +99,6 @@ class sparse_GP(GP):
# Compute dL_dKmm
tmp = tdot(self._LBi_Lmi_psi1V)
self.DBi_plus_BiPBi = backsub_both_sides(self.LB, self.D * np.eye(self.M) + tmp)
# tmp = -0.5 * self.DBi_plus_BiPBi / sf2
# tmp += -0.5 * self.B * sf2 * self.D
tmp = -0.5 * self.DBi_plus_BiPBi
tmp += -0.5 * self.B * self.D
tmp += self.D * np.eye(self.M)
@ -115,6 +108,7 @@ class sparse_GP(GP):
self.dL_dpsi0 = -0.5 * self.D * (self.likelihood.precision * np.ones([self.N, 1])).flatten()
self.dL_dpsi1 = np.dot(self.Cpsi1V, self.likelihood.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, :, :]
@ -141,7 +135,6 @@ class sparse_GP(GP):
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 * (self.psi0.sum() * self.likelihood.precision ** 2 - np.trace(self.A) * self.likelihood.precision)
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)))
@ -149,16 +142,12 @@ class sparse_GP(GP):
def log_likelihood(self):
""" Compute the (lower bound on the) log marginal likelihood """
# sf2 = self.scale_factor ** 2
if self.likelihood.is_heteroscedastic:
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.likelihood.V * self.likelihood.Y)
# B = -0.5 * self.D * (np.sum(self.likelihood.precision.flatten() * self.psi0) - np.trace(self.A) * sf2)
B = -0.5 * self.D * (np.sum(self.likelihood.precision.flatten() * self.psi0) - np.trace(self.A))
else:
A = -0.5 * self.N * self.D * (np.log(2.*np.pi) - np.log(self.likelihood.precision)) - 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)
B = -0.5 * self.D * (np.sum(self.likelihood.precision * self.psi0) - np.trace(self.A))
# C = -self.D * (np.sum(np.log(np.diag(self.LB))) + 0.5 * 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
@ -168,14 +157,6 @@ class sparse_GP(GP):
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()
# if self.auto_scale_factor:
# self.scale_factor = np.sqrt(self.psi2.sum(0).mean()*self.likelihood.precision)
# if self.auto_scale_factor:
# if self.likelihood.is_heteroscedastic:
# self.scale_factor = max(100,np.sqrt(self.psi2_beta_scaled.sum(0).mean()))
# else:
# self.scale_factor = np.sqrt(self.psi2.sum(0).mean()*self.likelihood.precision)
# self.scale_factor = 100.
self._computations()
def _get_params(self):
@ -188,7 +169,7 @@ class sparse_GP(GP):
"""
Approximates a non-gaussian likelihood using Expectation Propagation
For a Gaussian (or direct: TODO) likelihood, no iteration is required:
For a Gaussian likelihood, no iteration is required:
this function does nothing
"""
if self.has_uncertain_inputs: