diff --git a/GPy/core/gp.py b/GPy/core/gp.py index d8d1a87a..81985773 100644 --- a/GPy/core/gp.py +++ b/GPy/core/gp.py @@ -11,6 +11,7 @@ from parameterization import ObservableArray from .. import likelihoods from ..likelihoods.gaussian import Gaussian from ..inference.latent_function_inference import exact_gaussian_inference +from parameterization.variational import VariationalPosterior class GP(Model): """ @@ -30,10 +31,10 @@ class GP(Model): super(GP, self).__init__(name) assert X.ndim == 2 - if isinstance(X, ObservableArray): - self.X = self.X = X + if isinstance(X, ObservableArray) or isinstance(X, VariationalPosterior): + self.X = X else: self.X = ObservableArray(X) - + self.num_data, self.input_dim = self.X.shape assert Y.ndim == 2 diff --git a/GPy/core/parameterization/variational.py b/GPy/core/parameterization/variational.py index d1c0faf8..ef4d974d 100644 --- a/GPy/core/parameterization/variational.py +++ b/GPy/core/parameterization/variational.py @@ -35,6 +35,8 @@ class VariationalPosterior(Parameterized): def __init__(self, means=None, variances=None, name=None, **kw): super(VariationalPosterior, self).__init__(name=name, **kw) self.mean = Param("mean", means) + self.ndim = self.mean.ndim + self.shape = self.mean.shape self.variance = Param("variance", variances, Logexp()) self.add_parameters(self.mean, self.variance) self.num_data, self.input_dim = self.mean.shape diff --git a/GPy/core/sparse_gp.py b/GPy/core/sparse_gp.py index bb3116ba..a826cdf7 100644 --- a/GPy/core/sparse_gp.py +++ b/GPy/core/sparse_gp.py @@ -7,7 +7,7 @@ from gp import GP from parameterization.param import Param from ..inference.latent_function_inference import var_dtc from .. import likelihoods -from parameterization.variational import NormalPosterior +from parameterization.variational import VariationalPosterior class SparseGP(GP): """ @@ -32,7 +32,7 @@ class SparseGP(GP): """ - def __init__(self, X, Y, Z, kernel, likelihood, inference_method=None, X_variance=None, name='sparse gp'): + def __init__(self, X, Y, Z, kernel, likelihood, inference_method=None, name='sparse gp'): #pick a sensible inference method if inference_method is None: @@ -45,25 +45,24 @@ class SparseGP(GP): self.Z = Param('inducing inputs', Z) self.num_inducing = Z.shape[0] - - self.q = NormalPosterior(X, X_variance) - - GP.__init__(self, self.q.mean, Y, kernel, likelihood, inference_method=inference_method, name=name) + + GP.__init__(self, X, Y, kernel, likelihood, inference_method=inference_method, name=name) + self.add_parameter(self.Z, index=0) self.parameters_changed() def has_uncertain_inputs(self): - return self.q.has_uncertain_inputs() + if isinstance(self.X, VariationalPosterior): + return True + else: + return False def parameters_changed(self): - if self.has_uncertain_inputs(): - self.posterior, self._log_marginal_likelihood, self.grad_dict = self.inference_method.inference_latent(self.kern, self.q, self.Z, self.likelihood, self.Y) - else: - self.posterior, self._log_marginal_likelihood, self.grad_dict = self.inference_method.inference(self.kern, self.X, self.X_variance, self.Z, self.likelihood, self.Y) + self.posterior, self._log_marginal_likelihood, self.grad_dict = self.inference_method.inference(self.kern, self.X, self.Z, self.likelihood, self.Y) self.likelihood.update_gradients(self.grad_dict.pop('partial_for_likelihood')) - if self.has_uncertain_inputs(): - self.kern.update_gradients_variational(posterior_variational=self.q, Z=self.Z, **self.grad_dict) - self.Z.gradient = self.kern.gradients_Z_variational(posterior_variational=self.q, Z=self.Z, **self.grad_dict) + if isinstance(self.X, VariationalPosterior): + self.kern.update_gradients_variational(posterior_variational=self.X, Z=self.Z, **self.grad_dict) + self.Z.gradient = self.kern.gradients_Z_variational(posterior_variational=self.X, Z=self.Z, **self.grad_dict) else: self.kern.update_gradients_sparse(X=self.X, Z=self.Z, **self.grad_dict) self.Z.gradient = self.kern.gradients_Z_sparse(X=self.X, Z=self.Z, **self.grad_dict) diff --git a/GPy/inference/latent_function_inference/var_dtc.py b/GPy/inference/latent_function_inference/var_dtc.py index 349cd72d..ef4484bd 100644 --- a/GPy/inference/latent_function_inference/var_dtc.py +++ b/GPy/inference/latent_function_inference/var_dtc.py @@ -3,6 +3,7 @@ from posterior import Posterior from ...util.linalg import jitchol, backsub_both_sides, tdot, dtrtrs, dtrtri, dpotri, dpotrs, symmetrify +from ...core.parameterization.variational import VariationalPosterior import numpy as np from ...util.misc import param_to_array log_2_pi = np.log(2*np.pi) @@ -23,13 +24,13 @@ class VarDTC(object): from ...util.caching import Cacher self.get_trYYT = Cacher(self._get_trYYT, 1) self.get_YYTfactor = Cacher(self._get_YYTfactor, 1) - + def _get_trYYT(self, Y): return param_to_array(np.sum(np.square(Y))) def _get_YYTfactor(self, Y): """ - find a matrix L which satisfies LLT = YYT. + find a matrix L which satisfies LLT = YYT. Note that L may have fewer columns than Y. """ @@ -38,28 +39,26 @@ class VarDTC(object): return param_to_array(Y) else: return jitchol(tdot(Y)) - + def get_VVTfactor(self, Y, prec): return Y * prec # TODO chache this, and make it effective - - def inference(self, kern, X, X_variance, Z, likelihood, Y): - """Inference for normal sparseGP""" - uncertain_inputs = False - psi0, psi1, psi2 = _compute_psi(kern, X, X_variance, Z, uncertain_inputs) - return self._inference(kern, psi0, psi1, psi2, Z, likelihood, Y, uncertain_inputs) - def inference_latent(self, kern, posterior_variational, Z, likelihood, Y): - """Inference for GPLVM with uncertain inputs""" - uncertain_inputs = True - psi0, psi1, psi2 = _compute_psi_latent(kern, posterior_variational, Z) - return self._inference(kern, psi0, psi1, psi2, Z, likelihood, Y, uncertain_inputs) - - def _inference(self, kern, psi0, psi1, psi2, Z, likelihood, Y, uncertain_inputs): + def inference(self, kern, X, Z, likelihood, Y): + if isinstance(X, VariationalPosterior): + uncertain_inputs = True + psi0 = kern.psi0(Z, X) + psi1 = kern.psi1(Z, X) + psi2 = kern.psi2(Z, X) + else: + uncertain_inputs = False + psi0 = kern.Kdiag(X) + psi1 = kern.K(X, Z) + psi2 = None #see whether we're using variational uncertain inputs - + _, output_dim = Y.shape - + #see whether we've got a different noise variance for each datum beta = 1./np.squeeze(likelihood.variance) @@ -69,16 +68,16 @@ class VarDTC(object): VVT_factor = beta*Y #VVT_factor = beta*Y trYYT = self.get_trYYT(Y) - + # do the inference: het_noise = beta.size < 1 num_inducing = Z.shape[0] num_data = Y.shape[0] # kernel computations, using BGPLVM notation - Kmm = kern.K(Z) - + Kmm = kern.K(Z) + Lm = jitchol(Kmm) - + # The rather complex computations of A if uncertain_inputs: if het_noise: @@ -124,33 +123,33 @@ class VarDTC(object): dL_dKmm = backsub_both_sides(Lm, delit) # derivatives of L w.r.t. psi - dL_dpsi0, dL_dpsi1, dL_dpsi2 = _compute_dL_dpsi(num_inducing, num_data, output_dim, beta, Lm, - VVT_factor, Cpsi1Vf, DBi_plus_BiPBi, + dL_dpsi0, dL_dpsi1, dL_dpsi2 = _compute_dL_dpsi(num_inducing, num_data, output_dim, beta, Lm, + VVT_factor, Cpsi1Vf, DBi_plus_BiPBi, psi1, het_noise, uncertain_inputs) - + # log marginal likelihood - log_marginal = _compute_log_marginal_likelihood(likelihood, num_data, output_dim, beta, het_noise, + log_marginal = _compute_log_marginal_likelihood(likelihood, num_data, output_dim, beta, het_noise, psi0, A, LB, trYYT, data_fit) - + #put the gradients in the right places - partial_for_likelihood = _compute_partial_for_likelihood(likelihood, - het_noise, uncertain_inputs, LB, - _LBi_Lmi_psi1Vf, DBi_plus_BiPBi, Lm, A, - psi0, psi1, beta, + partial_for_likelihood = _compute_partial_for_likelihood(likelihood, + het_noise, uncertain_inputs, LB, + _LBi_Lmi_psi1Vf, DBi_plus_BiPBi, Lm, A, + psi0, psi1, beta, data_fit, num_data, output_dim, trYYT) - + #likelihood.update_gradients(partial_for_likelihood) if uncertain_inputs: - grad_dict = {'dL_dKmm': dL_dKmm, - 'dL_dpsi0':dL_dpsi0, - 'dL_dpsi1':dL_dpsi1, - 'dL_dpsi2':dL_dpsi2, + grad_dict = {'dL_dKmm': dL_dKmm, + 'dL_dpsi0':dL_dpsi0, + 'dL_dpsi1':dL_dpsi1, + 'dL_dpsi2':dL_dpsi2, 'partial_for_likelihood':partial_for_likelihood} else: - grad_dict = {'dL_dKmm': dL_dKmm, - 'dL_dKdiag':dL_dpsi0, - 'dL_dKnm':dL_dpsi1, + grad_dict = {'dL_dKmm': dL_dKmm, + 'dL_dKdiag':dL_dpsi0, + 'dL_dKnm':dL_dpsi1, 'partial_for_likelihood':partial_for_likelihood} #get sufficient things for posterior prediction @@ -181,7 +180,7 @@ class VarDTCMissingData(object): from ...util.caching import Cacher self._Y = Cacher(self._subarray_computations, 1) pass - + def _subarray_computations(self, Y): inan = np.isnan(Y) has_none = inan.any() @@ -202,19 +201,19 @@ class VarDTCMissingData(object): self._subarray_indices = [[slice(None),slice(None)]] return [Y], [(Y**2).sum()] - def inference(self, kern, X, X_variance, Z, likelihood, Y): - """Inference for normal sparseGP""" - uncertain_inputs = False - psi0, psi1, psi2 = _compute_psi(kern, X, X_variance, Z, uncertain_inputs) - return self._inference(kern, psi0, psi1, psi2, Z, likelihood, Y, uncertain_inputs) + def inference(self, kern, X, Z, likelihood, Y): + if isinstance(X, VariationalPosterior): + uncertain_inputs = True + psi0 = kern.psi0(Z, X) + psi1 = kern.psi1(Z, X) + psi2 = kern.psi2(Z, X) + else: + uncertain_inputs = False + psi0 = kern.Kdiag(X) + psi1 = kern.K(X, Z) + psi2 = None + - def inference_latent(self, kern, posterior_variational, Z, likelihood, Y): - """Inference for GPLVM with uncertain inputs""" - uncertain_inputs = True - psi0, psi1, psi2 = _compute_psi_latent(kern, posterior_variational, Z) - return self._inference(kern, psi0, psi1, psi2, Z, likelihood, Y, uncertain_inputs) - - def _inference(self, kern, psi0_all, psi1_all, psi2_all, Z, likelihood, Y, uncertain_inputs): Ys, traces = self._Y(Y) beta_all = 1./likelihood.variance het_noise = beta_all.size != 1 @@ -226,15 +225,15 @@ class VarDTCMissingData(object): dL_dpsi1_all = np.zeros((Y.shape[0], num_inducing)) if uncertain_inputs: dL_dpsi2_all = np.zeros((Y.shape[0], num_inducing, num_inducing)) - + partial_for_likelihood = 0 woodbury_vector = np.zeros((num_inducing, Y.shape[1])) woodbury_inv_all = np.zeros((num_inducing, num_inducing, Y.shape[1])) dL_dKmm = 0 log_marginal = 0 - + Kmm = kern.K(Z) - #factor Kmm + #factor Kmm Lm = jitchol(Kmm) if uncertain_inputs: LmInv = dtrtri(Lm) @@ -242,11 +241,11 @@ class VarDTCMissingData(object): full_VVT_factor = VVT_factor_all.shape[1] == Y.shape[1] if not full_VVT_factor: psi1V = np.dot(Y.T*beta_all, psi1_all).T - + for y, trYYT, [v, ind] in itertools.izip(Ys, traces, self._subarray_indices): if het_noise: beta = beta_all[ind] else: beta = beta_all[0] - + VVT_factor = (beta*y) VVT_factor_all[v, ind].flat = VVT_factor.flat output_dim = y.shape[1] @@ -256,7 +255,7 @@ class VarDTCMissingData(object): if uncertain_inputs: psi2 = psi2_all[v, :] else: psi2 = None num_data = psi1.shape[0] - + if uncertain_inputs: if het_noise: psi2_beta = psi2 * (beta.flatten().reshape(num_data, 1, 1)).sum(0) else: psi2_beta = psi2.sum(0) * beta @@ -270,13 +269,13 @@ class VarDTCMissingData(object): # factor B B = np.eye(num_inducing) + A LB = jitchol(B) - + psi1Vf = psi1.T.dot(VVT_factor) tmp, _ = dtrtrs(Lm, psi1Vf, lower=1, trans=0) _LBi_Lmi_psi1Vf, _ = dtrtrs(LB, tmp, lower=1, trans=0) tmp, _ = dtrtrs(LB, _LBi_Lmi_psi1Vf, lower=1, trans=1) Cpsi1Vf, _ = dtrtrs(Lm, tmp, lower=1, trans=1) - + # data fit and derivative of L w.r.t. Kmm delit = tdot(_LBi_Lmi_psi1Vf) data_fit = np.trace(delit) @@ -287,34 +286,34 @@ class VarDTCMissingData(object): dL_dKmm += backsub_both_sides(Lm, delit) # derivatives of L w.r.t. psi - dL_dpsi0, dL_dpsi1, dL_dpsi2 = _compute_dL_dpsi(num_inducing, num_data, output_dim, beta, Lm, - VVT_factor, Cpsi1Vf, DBi_plus_BiPBi, + dL_dpsi0, dL_dpsi1, dL_dpsi2 = _compute_dL_dpsi(num_inducing, num_data, output_dim, beta, Lm, + VVT_factor, Cpsi1Vf, DBi_plus_BiPBi, psi1, het_noise, uncertain_inputs) - + #import ipdb;ipdb.set_trace() dL_dpsi0_all[v] += dL_dpsi0 dL_dpsi1_all[v, :] += dL_dpsi1 if uncertain_inputs: dL_dpsi2_all[v, :] += dL_dpsi2 - + # log marginal likelihood - log_marginal += _compute_log_marginal_likelihood(likelihood, num_data, output_dim, beta, het_noise, + log_marginal += _compute_log_marginal_likelihood(likelihood, num_data, output_dim, beta, het_noise, psi0, A, LB, trYYT, data_fit) #put the gradients in the right places - partial_for_likelihood += _compute_partial_for_likelihood(likelihood, - het_noise, uncertain_inputs, LB, - _LBi_Lmi_psi1Vf, DBi_plus_BiPBi, Lm, A, - psi0, psi1, beta, + partial_for_likelihood += _compute_partial_for_likelihood(likelihood, + het_noise, uncertain_inputs, LB, + _LBi_Lmi_psi1Vf, DBi_plus_BiPBi, Lm, A, + psi0, psi1, beta, data_fit, num_data, output_dim, trYYT) - + if full_VVT_factor: woodbury_vector[:, ind] = Cpsi1Vf else: print 'foobar' tmp, _ = dtrtrs(Lm, psi1V, lower=1, trans=0) tmp, _ = dpotrs(LB, tmp, lower=1) woodbury_vector[:, ind] = dtrtrs(Lm, tmp, lower=1, trans=1)[0] - + #import ipdb;ipdb.set_trace() Bi, _ = dpotri(LB, lower=1) symmetrify(Bi) @@ -325,15 +324,15 @@ class VarDTCMissingData(object): # gradients: if uncertain_inputs: - grad_dict = {'dL_dKmm': dL_dKmm, - 'dL_dpsi0':dL_dpsi0_all, - 'dL_dpsi1':dL_dpsi1_all, - 'dL_dpsi2':dL_dpsi2_all, + grad_dict = {'dL_dKmm': dL_dKmm, + 'dL_dpsi0':dL_dpsi0_all, + 'dL_dpsi1':dL_dpsi1_all, + 'dL_dpsi2':dL_dpsi2_all, 'partial_for_likelihood':partial_for_likelihood} else: - grad_dict = {'dL_dKmm': dL_dKmm, - 'dL_dKdiag':dL_dpsi0_all, - 'dL_dKnm':dL_dpsi1_all, + grad_dict = {'dL_dKmm': dL_dKmm, + 'dL_dKdiag':dL_dpsi0_all, + 'dL_dKnm':dL_dpsi1_all, 'partial_for_likelihood':partial_for_likelihood} #get sufficient things for posterior prediction @@ -350,26 +349,13 @@ class VarDTCMissingData(object): #Bi = -dpotri(LB_all, lower=1)[0] #from ...util import diag #diag.add(Bi, 1) - + #woodbury_inv = backsub_both_sides(Lm, Bi) - + post = Posterior(woodbury_inv=woodbury_inv_all, woodbury_vector=woodbury_vector, K=Kmm, mean=None, cov=None, K_chol=Lm) return post, log_marginal, grad_dict - -def _compute_psi(kern, X, X_variance, Z): - psi0 = kern.Kdiag(X) - psi1 = kern.K(X, Z) - psi2 = None - return psi0, psi1, psi2 - -def _compute_psi_latent(kern, posterior_variational, Z): - psi0 = kern.psi0(Z, posterior_variational) - psi1 = kern.psi1(Z, posterior_variational) - psi2 = kern.psi2(Z, posterior_variational) - return psi0, psi1, psi2 - def _compute_dL_dpsi(num_inducing, num_data, output_dim, beta, Lm, VVT_factor, Cpsi1Vf, DBi_plus_BiPBi, psi1, het_noise, uncertain_inputs): dL_dpsi0 = -0.5 * output_dim * (beta * np.ones([num_data, 1])).flatten() dL_dpsi1 = np.dot(VVT_factor, Cpsi1Vf.T) diff --git a/GPy/models/bayesian_gplvm.py b/GPy/models/bayesian_gplvm.py index 74b8abe0..50fc2810 100644 --- a/GPy/models/bayesian_gplvm.py +++ b/GPy/models/bayesian_gplvm.py @@ -32,21 +32,23 @@ class BayesianGPLVM(SparseGP): if X_variance is None: X_variance = np.random.uniform(0,.1,X.shape) + if Z is None: Z = np.random.permutation(X.copy())[:num_inducing] assert Z.shape[1] == X.shape[1] if kernel is None: kernel = kern.RBF(input_dim) # + kern.white(input_dim) - + if likelihood is None: likelihood = Gaussian() - self.q = NormalPosterior(X, X_variance) + + self.variational_prior = NormalPrior() - - SparseGP.__init__(self, X, Y, Z, kernel, likelihood, inference_method, X_variance, name, **kwargs) - self.add_parameter(self.q, index=0) - #self.ensure_default_constraints() + X = NormalPosterior(X, X_variance) + + SparseGP.__init__(self, X, Y, Z, kernel, likelihood, inference_method, name, **kwargs) + self.add_parameter(self.X, index=0) def _getstate(self): """ @@ -62,16 +64,14 @@ class BayesianGPLVM(SparseGP): def parameters_changed(self): super(BayesianGPLVM, self).parameters_changed() - self._log_marginal_likelihood -= self.variational_prior.KL_divergence(self.q) - - # TODO: This has to go into kern - # maybe a update_gradients_q_variational? - self.q.mean.gradient, self.q.variance.gradient = self.kern.gradients_q_variational(posterior_variational=self.q, Z=self.Z, **self.grad_dict) - + self._log_marginal_likelihood -= self.variational_prior.KL_divergence(self.X) + + self.X.mean.gradient, self.X.variance.gradient = self.kern.gradients_q_variational(posterior_variational=self.X, Z=self.Z, **self.grad_dict) + # update for the KL divergence - self.variational_prior.update_gradients_KL(self.q) - - + self.variational_prior.update_gradients_KL(self.X) + + def plot_latent(self, plot_inducing=True, *args, **kwargs): """ See GPy.plotting.matplot_dep.dim_reduction_plots.plot_latent @@ -150,14 +150,14 @@ class BayesianGPLVM(SparseGP): return dim_reduction_plots.plot_steepest_gradient_map(self,*args,**kwargs) class BayesianGPLVMWithMissingData(BayesianGPLVM): - def __init__(self, Y, input_dim, X=None, X_variance=None, init='PCA', num_inducing=10, + def __init__(self, Y, input_dim, X=None, X_variance=None, init='PCA', num_inducing=10, Z=None, kernel=None, inference_method=None, likelihood=None, name='bayesian gplvm', **kwargs): from ..util.subarray_and_sorting import common_subarrays self.subarrays = common_subarrays(Y) import ipdb;ipdb.set_trace() BayesianGPLVM.__init__(self, Y, input_dim, X=X, X_variance=X_variance, init=init, num_inducing=num_inducing, Z=Z, kernel=kernel, inference_method=inference_method, likelihood=likelihood, name=name, **kwargs) - - + + def parameters_changed(self): super(BayesianGPLVM, self).parameters_changed() self._log_marginal_likelihood -= self.KL_divergence() @@ -165,12 +165,12 @@ class BayesianGPLVMWithMissingData(BayesianGPLVM): dL_dmu, dL_dS = self.dL_dmuS() # dL: - self.q.mean.gradient = dL_dmu - self.q.variance.gradient = dL_dS + self.X.mean.gradient = dL_dmu + self.X.variance.gradient = dL_dS # dKL: - self.q.mean.gradient -= self.X - self.q.variance.gradient -= (1. - (1. / (self.X_variance))) * 0.5 + self.X.mean.gradient -= self.X.mean + self.X.variance.gradient -= (1. - (1. / (self.X.variance))) * 0.5 if __name__ == '__main__': import numpy as np @@ -178,7 +178,7 @@ if __name__ == '__main__': W = np.linspace(0,1,10)[None,:] Y = (X*W).sum(1) missing = np.random.binomial(1,.1,size=Y.shape) - + pass def latent_cost_and_grad(mu_S, kern, Z, dL_dpsi0, dL_dpsi1, dL_dpsi2):