From 78b50db138beff11ea528e6155d2a1bd28554d2c Mon Sep 17 00:00:00 2001 From: Zhenwen Dai Date: Mon, 3 Nov 2014 13:38:28 +0000 Subject: [PATCH] add new inference X --- GPy/core/parameterization/variational.py | 6 + .../latent_function_inference/inferenceX.py | 127 ++++++++++++++++++ GPy/models/bayesian_gplvm.py | 16 +++ 3 files changed, 149 insertions(+) create mode 100644 GPy/inference/latent_function_inference/inferenceX.py diff --git a/GPy/core/parameterization/variational.py b/GPy/core/parameterization/variational.py index 251ec7db..e2a24008 100644 --- a/GPy/core/parameterization/variational.py +++ b/GPy/core/parameterization/variational.py @@ -94,6 +94,9 @@ class VariationalPosterior(Parameterized): if self.has_uncertain_inputs(): assert self.variance.shape == self.mean.shape, "need one variance per sample and dimenion" + def set_gradients(self, grad): + self.mean.gradient, self.variance.gradient = grad + def _raveled_index(self): index = np.empty(dtype=int, shape=0) size = 0 @@ -158,6 +161,9 @@ class SpikeAndSlabPosterior(VariationalPosterior): self.gamma = Param("binary_prob",binary_prob, Logistic(1e-10,1.-1e-10)) self.link_parameter(self.gamma) + def set_gradients(self, grad): + self.mean.gradient, self.variance.gradient, self.gamma.gradient = grad + def __getitem__(self, s): if isinstance(s, (int, slice, tuple, list, np.ndarray)): import copy diff --git a/GPy/inference/latent_function_inference/inferenceX.py b/GPy/inference/latent_function_inference/inferenceX.py new file mode 100644 index 00000000..9828c4a1 --- /dev/null +++ b/GPy/inference/latent_function_inference/inferenceX.py @@ -0,0 +1,127 @@ +""" +""" +import numpy as np +from ...core import Model +from ...core.parameterization import variational + +def infer_newX(model, Y_new, optimize=True, init='L2'): + """ + Infer the distribution of X for the new observed data *Y_new*. + + :param model: the GPy model used in inference + :type model: GPy.core.Model + :param Y_new: the new observed data for inference + :type Y_new: numpy.ndarray + :param optimize: whether to optimize the location of new X (True by default) + :type optimize: boolean + :return: a tuple containing the estimated posterior distribution of X and the model that optimize X + :rtype: (GPy.core.parameterization.variational.VariationalPosterior, GPy.core.Model) + """ + infr_m = InferenceX(model, Y_new, init=init) + + if optimize: + infr_m.optimize() + + return infr_m.X, infr_m + +class InferenceX(Model): + """ + The class for inference of new X with given new Y. (do_test_latent) + + :param model: the GPy model used in inference + :type model: GPy.core.Model + :param Y: the new observed data for inference + :type Y: numpy.ndarray + """ + def __init__(self, model, Y, name='inferenceX', init='L2'): + if np.isnan(Y).any(): + assert Y.shape[0]==1, "The current implementation of inference X only support one data point at a time with missing data!" + self.missing_data = True + self.valid_dim = np.logical_not(np.isnan(Y[0])) + else: + self.missing_data = False + super(InferenceX, self).__init__(name) + self.likelihood = model.likelihood.copy() + self.kern = model.kern.copy() + if model.kern.useGPU: + from ...models import SSGPLVM + if isinstance(model, SSGPLVM): + self.kern.GPU_SSRBF(True) + else: + self.kern.GPU(True) + from copy import deepcopy + self.posterior = deepcopy(model.posterior) + self.variational_prior = model.variational_prior.copy() + self.Z = model.Z.copy() + self.Y = Y + self.X = self._init_X(model, Y, init=init) + self.compute_dL() + + self.link_parameter(self.X) + + def _init_X(self, model, Y_new, init='L2'): + # Initialize the new X by finding the nearest point in Y space. + + Y = model.Y + if self.missing_data: + Y = Y[:,self.valid_dim] + Y_new = Y_new[:,self.valid_dim] + dist = -2.*Y_new.dot(Y.T) + np.square(Y_new).sum(axis=1)[:,None]+ np.square(Y).sum(axis=1)[None,:] + else: + if init=='L2': + dist = -2.*Y_new.dot(Y.T) + np.square(Y_new).sum(axis=1)[:,None]+ np.square(Y).sum(axis=1)[None,:] + elif init=='NCC': + dist = Y_new.dot(Y.T) + idx = dist.argmin(axis=1) + + from ...models import SSGPLVM + from ...util.misc import param_to_array + if isinstance(model, SSGPLVM): + X = variational.SpikeAndSlabPosterior(param_to_array(model.X.mean[idx]), param_to_array(model.X.variance[idx]), param_to_array(model.X.gamma[idx])) + if model.group_spike: + X.gamma.fix() + else: + X = variational.NormalPosterior(param_to_array(model.X.mean[idx]), param_to_array(model.X.variance[idx])) + + return X + + def compute_dL(self): + # Common computation + beta = 1./np.fmax(self.likelihood.variance, 1e-6) + output_dim = self.Y.shape[-1] + wv = self.posterior.woodbury_vector + if self.missing_data: + wv = wv[:,self.valid_dim] + output_dim = self.valid_dim.sum() + self.dL_dpsi2 = beta*(output_dim*self.posterior.woodbury_inv - np.einsum('md,od->mo',wv, wv))/2. + self.dL_dpsi1 = beta*np.dot(self.Y[:,self.valid_dim], wv.T) + self.dL_dpsi0 = -output_dim * beta/2.* np.ones(self.Y.shape[0]) + else: + self.dL_dpsi2 = beta*(output_dim*self.posterior.woodbury_inv - np.einsum('md,od->mo',wv, wv))/2. + self.dL_dpsi1 = beta*np.dot(self.Y, wv.T) + self.dL_dpsi0 = -output_dim * beta/2.* np.ones(self.Y.shape[0]) + + def parameters_changed(self): + psi0 = self.kern.psi0(self.Z, self.X) + psi1 = self.kern.psi1(self.Z, self.X) + psi2 = self.kern.psi2(self.Z, self.X) + + self._log_marginal_likelihood = (self.dL_dpsi2*psi2).sum()+(self.dL_dpsi1*psi1).sum()+(self.dL_dpsi0*psi0).sum() + X_grad = self.kern.gradients_qX_expectations(variational_posterior=self.X, Z=self.Z, dL_dpsi0=self.dL_dpsi0, dL_dpsi1=self.dL_dpsi1, dL_dpsi2=self.dL_dpsi2) + self.X.set_gradients(X_grad) + + from ...core.parameterization.variational import SpikeAndSlabPrior + if isinstance(self.variational_prior, SpikeAndSlabPrior): + # Update Log-likelihood + KL_div = self.variational_prior.KL_divergence(self.X, N=self.Y.shape[0]) + # update for the KL divergence + self.variational_prior.update_gradients_KL(self.X, N=self.Y.shape[0]) + else: + # Update Log-likelihood + KL_div = self.variational_prior.KL_divergence(self.X) + # update for the KL divergence + self.variational_prior.update_gradients_KL(self.X) + + def log_likelihood(self): + return self._log_marginal_likelihood + diff --git a/GPy/models/bayesian_gplvm.py b/GPy/models/bayesian_gplvm.py index c4d0b6e9..b2eed718 100644 --- a/GPy/models/bayesian_gplvm.py +++ b/GPy/models/bayesian_gplvm.py @@ -175,6 +175,22 @@ class BayesianGPLVM(SparseGP_MPI): resolution, ax, marker, s, fignum, plot_inducing, legend, plot_limits, aspect, updates, predict_kwargs, imshow_kwargs) + + def infer_newX(self, Y_new, optimize=True, ): + """ + Infer the distribution of X for the new observed data *Y_new*. + + :param model: the GPy model used in inference + :type model: GPy.core.Model + :param Y_new: the new observed data for inference + :type Y_new: numpy.ndarray + :param optimize: whether to optimize the location of new X (True by default) + :type optimize: boolean + :return: a tuple containing the estimated posterior distribution of X and the model that optimize X + :rtype: (GPy.core.parameterization.variational.VariationalPosterior, GPy.core.Model) + """ + from ..inference.latent_function_inference.inferenceX import infer_newX + return infer_newX(self, Y_new, optimize=optimize) def do_test_latents(self, Y): """