From 6b3888f1631bd94f7cd7b8010946d42a960cef76 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Mon, 3 Nov 2014 11:16:34 +0000 Subject: [PATCH 1/5] [VarDTC] reverted SparseGP to previous state, updated BGPLVM accordingly --- GPy/core/sparse_gp.py | 282 +++----------------- GPy/core/sparse_gp_mpi.py | 6 +- GPy/models/bayesian_gplvm_minibatch.py | 267 +++++++++++++++++++ GPy/models/sparse_gp_minibatch.py | 347 +++++++++++++++++++++++++ GPy/testing/model_tests.py | 6 +- 5 files changed, 655 insertions(+), 253 deletions(-) create mode 100644 GPy/models/bayesian_gplvm_minibatch.py create mode 100644 GPy/models/sparse_gp_minibatch.py diff --git a/GPy/core/sparse_gp.py b/GPy/core/sparse_gp.py index 7f88503a..73c80c76 100644 --- a/GPy/core/sparse_gp.py +++ b/GPy/core/sparse_gp.py @@ -40,8 +40,7 @@ class SparseGP(GP): """ def __init__(self, X, Y, Z, kernel, likelihood, inference_method=None, - name='sparse gp', Y_metadata=None, normalizer=False, - missing_data=False, stochastic=False, batchsize=1): + name='sparse gp', Y_metadata=None, normalizer=False): #pick a sensible inference method if inference_method is None: if isinstance(likelihood, likelihoods.Gaussian): @@ -55,260 +54,51 @@ class SparseGP(GP): self.num_inducing = Z.shape[0] GP.__init__(self, X, Y, kernel, likelihood, inference_method=inference_method, name=name, Y_metadata=Y_metadata, normalizer=normalizer) - self.missing_data = missing_data - - if stochastic and missing_data: - self.missing_data = True - self.ninan = ~np.isnan(Y) - self.stochastics = SparseGPStochastics(self, batchsize) - elif stochastic and not missing_data: - self.missing_data = False - self.stochastics = SparseGPStochastics(self, batchsize) - elif missing_data: - self.missing_data = True - self.ninan = ~np.isnan(Y) - self.stochastics = SparseGPMissing(self) - else: - self.stochastics = False logger.info("Adding Z as parameter") self.link_parameter(self.Z, index=0) - if self.missing_data: - self.Ylist = [] - overall = self.Y_normalized.shape[1] - m_f = lambda i: "Precomputing Y for missing data: {: >7.2%}".format(float(i+1)/overall) - message = m_f(-1) - print message, - for d in xrange(overall): - self.Ylist.append(self.Y_normalized[self.ninan[:, d], d][:, None]) - print ' '*(len(message)+1) + '\r', - message = m_f(d) - print message, - print '' - self.posterior = None def has_uncertain_inputs(self): return isinstance(self.X, VariationalPosterior) - def _inner_parameters_changed(self, kern, X, Z, likelihood, Y, Y_metadata, Lm=None, dL_dKmm=None, subset_indices=None): - """ - This is the standard part, which usually belongs in parameters_changed. - - For automatic handling of subsampling (such as missing_data, stochastics etc.), we need to put this into an inner - loop, in order to ensure a different handling of gradients etc of different - subsets of data. - - The dict in current_values will be passed aroung as current_values for - the rest of the algorithm, so this is the place to store current values, - such as subsets etc, if necessary. - - If Lm and dL_dKmm can be precomputed (or only need to be computed once) - pass them in here, so they will be passed to the inference_method. - - subset_indices is a dictionary of indices. you can put the indices however you - like them into this dictionary for inner use of the indices inside the - algorithm. - """ - try: - posterior, log_marginal_likelihood, grad_dict = self.inference_method.inference(kern, X, Z, likelihood, Y, Y_metadata, Lm=Lm, dL_dKmm=None) - except: - posterior, log_marginal_likelihood, grad_dict = self.inference_method.inference(kern, X, Z, likelihood, Y, Y_metadata) - current_values = {} - likelihood.update_gradients(grad_dict['dL_dthetaL']) - current_values['likgrad'] = likelihood.gradient.copy() - if subset_indices is None: - subset_indices = {} - if isinstance(X, VariationalPosterior): - #gradients wrt kernel - dL_dKmm = grad_dict['dL_dKmm'] - kern.update_gradients_full(dL_dKmm, Z, None) - current_values['kerngrad'] = kern.gradient.copy() - kern.update_gradients_expectations(variational_posterior=X, - Z=Z, - dL_dpsi0=grad_dict['dL_dpsi0'], - dL_dpsi1=grad_dict['dL_dpsi1'], - dL_dpsi2=grad_dict['dL_dpsi2']) - current_values['kerngrad'] += kern.gradient - - #gradients wrt Z - current_values['Zgrad'] = kern.gradients_X(dL_dKmm, Z) - current_values['Zgrad'] += kern.gradients_Z_expectations( - grad_dict['dL_dpsi0'], - grad_dict['dL_dpsi1'], - grad_dict['dL_dpsi2'], - Z=Z, - variational_posterior=X) - else: - #gradients wrt kernel - kern.update_gradients_diag(grad_dict['dL_dKdiag'], X) - current_values['kerngrad'] = kern.gradient.copy() - kern.update_gradients_full(grad_dict['dL_dKnm'], X, Z) - current_values['kerngrad'] += kern.gradient - kern.update_gradients_full(grad_dict['dL_dKmm'], Z, None) - current_values['kerngrad'] += kern.gradient - #gradients wrt Z - current_values['Zgrad'] = kern.gradients_X(grad_dict['dL_dKmm'], Z) - current_values['Zgrad'] += kern.gradients_X(grad_dict['dL_dKnm'].T, Z, X) - return posterior, log_marginal_likelihood, grad_dict, current_values, subset_indices - - def _inner_take_over_or_update(self, full_values=None, current_values=None, value_indices=None): - """ - This is for automatic updates of values in the inner loop of missing - data handling. Both arguments are dictionaries and the values in - full_values will be updated by the current_gradients. - - If a key from current_values does not exist in full_values, it will be - initialized to the value in current_values. - - If there is indices needed for the update, value_indices can be used for - that. If value_indices has the same key, as current_values, the update - in full_values will be indexed by the indices in value_indices. - - grads: - dictionary of standing gradients (you will have to carefully make sure, that - the ordering is right!). The values in here will be updated such that - full_values[key] += current_values[key] forall key in full_gradients.keys() - - gradients: - dictionary of gradients in the current set of parameters. - - value_indices: - dictionary holding indices for the update in full_values. - if the key exists the update rule is:def df(x): - full_values[key][value_indices[key]] += current_values[key] - """ - for key in current_values.keys(): - if value_indices is not None and value_indices.has_key(key): - index = value_indices[key] - else: - index = slice(None) - if full_values.has_key(key): - full_values[key][index] += current_values[key] - else: - full_values[key] = current_values[key] - - def _inner_values_update(self, current_values): - """ - This exists if there is more to do with the current values. - It will be called allways in the inner loop, so that - you can do additional inner updates for the inside of the missing data - loop etc. This can also be used for stochastic updates, when only working on - one dimension of the output. - """ - pass - - def _outer_values_update(self, full_values): - """ - Here you put the values, which were collected before in the right places. - E.g. set the gradients of parameters, etc. - """ - self.likelihood.gradient = full_values['likgrad'] - self.kern.gradient = full_values['kerngrad'] - self.Z.gradient = full_values['Zgrad'] - - def _outer_init_full_values(self): - """ - If full_values has indices in values_indices, we might want to initialize - the full_values differently, so that subsetting is possible. - - Here you can initialize the full_values for the values needed. - - Keep in mind, that if a key does not exist in full_values when updating - values, it will be set (so e.g. for Z there is no need to initialize Zgrad, - as there is no subsetting needed. For X in BGPLVM on the other hand we probably need - to initialize the gradients for the mean and the variance in order to - have the full gradient for indexing) - """ - return {} - - def _outer_loop_for_missing_data(self): - Lm = None - dL_dKmm = None - - self._log_marginal_likelihood = 0 - self.full_values = self._outer_init_full_values() - - if self.posterior is None: - woodbury_inv = np.zeros((self.num_inducing, self.num_inducing, self.output_dim)) - woodbury_vector = np.zeros((self.num_inducing, self.output_dim)) - else: - woodbury_inv = self.posterior._woodbury_inv - woodbury_vector = self.posterior._woodbury_vector - - if not self.stochastics: - m_f = lambda i: "Inference with missing_data: {: >7.2%}".format(float(i+1)/self.output_dim) - message = m_f(-1) - print message, - - for d in self.stochastics.d: - ninan = self.ninan[:, d] - - if not self.stochastics: - print ' '*(len(message)) + '\r', - message = m_f(d) - print message, - - posterior, log_marginal_likelihood, \ - grad_dict, current_values, value_indices = self._inner_parameters_changed( - self.kern, self.X[ninan], - self.Z, self.likelihood, - self.Ylist[d], self.Y_metadata, - Lm, dL_dKmm, - subset_indices=dict(outputs=d, samples=ninan)) - - self._inner_take_over_or_update(self.full_values, current_values, value_indices) - self._inner_values_update(current_values) - - Lm = posterior.K_chol - dL_dKmm = grad_dict['dL_dKmm'] - woodbury_inv[:, :, d] = posterior.woodbury_inv - woodbury_vector[:, d:d+1] = posterior.woodbury_vector - self._log_marginal_likelihood += log_marginal_likelihood - if not self.stochastics: - print '' - - if self.posterior is None: - self.posterior = Posterior(woodbury_inv=woodbury_inv, woodbury_vector=woodbury_vector, - K=posterior._K, mean=None, cov=None, K_chol=posterior.K_chol) - self._outer_values_update(self.full_values) - - def _outer_loop_without_missing_data(self): - self._log_marginal_likelihood = 0 - - if self.posterior is None: - woodbury_inv = np.zeros((self.num_inducing, self.num_inducing, self.output_dim)) - woodbury_vector = np.zeros((self.num_inducing, self.output_dim)) - else: - woodbury_inv = self.posterior._woodbury_inv - woodbury_vector = self.posterior._woodbury_vector - - d = self.stochastics.d - posterior, log_marginal_likelihood, \ - grad_dict, self.full_values, _ = self._inner_parameters_changed( - self.kern, self.X, - self.Z, self.likelihood, - self.Y_normalized[:, d], self.Y_metadata) - self.grad_dict = grad_dict - - self._log_marginal_likelihood += log_marginal_likelihood - - self._outer_values_update(self.full_values) - - woodbury_inv[:, :, d] = posterior.woodbury_inv[:, :, None] - woodbury_vector[:, d] = posterior.woodbury_vector - if self.posterior is None: - self.posterior = Posterior(woodbury_inv=woodbury_inv, woodbury_vector=woodbury_vector, - K=posterior._K, mean=None, cov=None, K_chol=posterior.K_chol) - def parameters_changed(self): - if self.missing_data: - self._outer_loop_for_missing_data() - elif self.stochastics: - self._outer_loop_without_missing_data() + self.posterior, self._log_marginal_likelihood, self.grad_dict = self.inference_method.inference(self.kern, self.X, self.Z, self.likelihood, self.Y, self.Y_metadata) + + self.likelihood.update_gradients(self.grad_dict['dL_dthetaL']) + + if isinstance(self.X, VariationalPosterior): + #gradients wrt kernel + dL_dKmm = self.grad_dict['dL_dKmm'] + self.kern.update_gradients_full(dL_dKmm, self.Z, None) + kerngrad = self.kern.gradient.copy() + self.kern.update_gradients_expectations(variational_posterior=self.X, + Z=self.Z, + dL_dpsi0=self.grad_dict['dL_dpsi0'], + dL_dpsi1=self.grad_dict['dL_dpsi1'], + dL_dpsi2=self.grad_dict['dL_dpsi2']) + self.kern.gradient += kerngrad + + #gradients wrt Z + self.Z.gradient = self.kern.gradients_X(dL_dKmm, self.Z) + self.Z.gradient += self.kern.gradients_Z_expectations( + self.grad_dict['dL_dpsi0'], + self.grad_dict['dL_dpsi1'], + self.grad_dict['dL_dpsi2'], + Z=self.Z, + variational_posterior=self.X) else: - self.posterior, self._log_marginal_likelihood, self.grad_dict, self.full_values, _ = self._inner_parameters_changed(self.kern, self.X, self.Z, self.likelihood, self.Y_normalized, self.Y_metadata) - self._outer_values_update(self.full_values) + #gradients wrt kernel + self.kern.update_gradients_diag(self.grad_dict['dL_dKdiag'], self.X) + kerngrad = self.kern.gradient.copy() + self.kern.update_gradients_full(self.grad_dict['dL_dKnm'], self.X, self.Z) + kerngrad += self.kern.gradient + self.kern.update_gradients_full(self.grad_dict['dL_dKmm'], self.Z, None) + self.kern.gradient += kerngrad + #gradients wrt Z + self.Z.gradient = self.kern.gradients_X(self.grad_dict['dL_dKmm'], self.Z) + self.Z.gradient += self.kern.gradients_X(self.grad_dict['dL_dKnm'].T, self.Z, self.X) + def _raw_predict(self, Xnew, full_cov=False, kern=None): """ diff --git a/GPy/core/sparse_gp_mpi.py b/GPy/core/sparse_gp_mpi.py index 52f75d45..733d0e1c 100644 --- a/GPy/core/sparse_gp_mpi.py +++ b/GPy/core/sparse_gp_mpi.py @@ -34,8 +34,7 @@ class SparseGP_MPI(SparseGP): """ - def __init__(self, X, Y, Z, kernel, likelihood, variational_prior=None, inference_method=None, name='sparse gp mpi', Y_metadata=None, mpi_comm=None, normalizer=False, - missing_data=False, stochastic=False, batchsize=1): + def __init__(self, X, Y, Z, kernel, likelihood, variational_prior=None, inference_method=None, name='sparse gp mpi', Y_metadata=None, mpi_comm=None, normalizer=False): self._IN_OPTIMIZATION_ = False if mpi_comm != None: if inference_method is None: @@ -43,8 +42,7 @@ class SparseGP_MPI(SparseGP): else: assert isinstance(inference_method, VarDTC_minibatch), 'inference_method has to support MPI!' - super(SparseGP_MPI, self).__init__(X, Y, Z, kernel, likelihood, inference_method=inference_method, name=name, Y_metadata=Y_metadata, normalizer=normalizer, - missing_data=missing_data, stochastic=stochastic, batchsize=batchsize) + super(SparseGP_MPI, self).__init__(X, Y, Z, kernel, likelihood, inference_method=inference_method, name=name, Y_metadata=Y_metadata, normalizer=normalizer) self.update_model(False) self.link_parameter(self.X, index=0) if variational_prior is not None: diff --git a/GPy/models/bayesian_gplvm_minibatch.py b/GPy/models/bayesian_gplvm_minibatch.py new file mode 100644 index 00000000..56b46244 --- /dev/null +++ b/GPy/models/bayesian_gplvm_minibatch.py @@ -0,0 +1,267 @@ +# Copyright (c) 2012, GPy authors (see AUTHORS.txt). +# Licensed under the BSD 3-clause license (see LICENSE.txt) + +import numpy as np +from .. import kern +from ..likelihoods import Gaussian +from ..core.parameterization.variational import NormalPosterior, NormalPrior +from ..inference.latent_function_inference.var_dtc_parallel import VarDTC_minibatch +from ..inference.latent_function_inference.var_dtc_gpu import VarDTC_GPU +import logging +from GPy.models.sparse_gp_minibatch import SparseGPMiniBatch + +class BayesianGPLVMMiniBatch(SparseGPMiniBatch): + """ + Bayesian Gaussian Process Latent Variable Model + + :param Y: observed data (np.ndarray) or GPy.likelihood + :type Y: np.ndarray| GPy.likelihood instance + :param input_dim: latent dimensionality + :type input_dim: int + :param init: initialisation method for the latent space + :type init: 'PCA'|'random' + + """ + 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', normalizer=None, + missing_data=False, stochastic=False, batchsize=1): + self.__IN_OPTIMIZATION__ = False + + self.logger = logging.getLogger(self.__class__.__name__) + if X is None: + from ..util.initialization import initialize_latent + self.logger.info("initializing latent space X with method {}".format(init)) + X, fracs = initialize_latent(init, input_dim, Y) + else: + fracs = np.ones(input_dim) + + self.init = init + + if X_variance is None: + self.logger.info("initializing latent space variance ~ uniform(0,.1)") + X_variance = np.random.uniform(0,.1,X.shape) + + if Z is None: + self.logger.info("initializing inducing inputs") + Z = np.random.permutation(X.copy())[:num_inducing] + assert Z.shape[1] == X.shape[1] + + if kernel is None: + self.logger.info("initializing kernel RBF") + kernel = kern.RBF(input_dim, lengthscale=1./fracs, ARD=True) #+ kern.Bias(input_dim) + kern.White(input_dim) + + if likelihood is None: + likelihood = Gaussian() + + self.variational_prior = NormalPrior() + X = NormalPosterior(X, X_variance) + + if inference_method is None: + from ..inference.latent_function_inference.var_dtc import VarDTC + self.logger.debug("creating inference_method var_dtc") + inference_method = VarDTC(limit=1 if not missing_data else Y.shape[1]) + + if kernel.useGPU and isinstance(inference_method, VarDTC_GPU): + kernel.psicomp.GPU_direct = True + + super(BayesianGPLVMMiniBatch,self).__init__(X, Y, Z, kernel, likelihood=likelihood, + name=name, inference_method=inference_method, + normalizer=normalizer, + missing_data=missing_data, stochastic=stochastic, + batchsize=batchsize) + + def set_X_gradients(self, X, X_grad): + """Set the gradients of the posterior distribution of X in its specific form.""" + X.mean.gradient, X.variance.gradient = X_grad + + def get_X_gradients(self, X): + """Get the gradients of the posterior distribution of X in its specific form.""" + return X.mean.gradient, X.variance.gradient + + def _inner_parameters_changed(self, kern, X, Z, likelihood, Y, Y_metadata, Lm=None, dL_dKmm=None, subset_indices=None): + posterior, log_marginal_likelihood, grad_dict, current_values, value_indices = super(BayesianGPLVMMiniBatch, self)._inner_parameters_changed(kern, X, Z, likelihood, Y, Y_metadata, Lm=Lm, dL_dKmm=dL_dKmm, subset_indices=subset_indices) + + kl_fctr = 1. + if self.missing_data: + d = self.output_dim + log_marginal_likelihood -= kl_fctr*self.variational_prior.KL_divergence(X)/d + else: + log_marginal_likelihood -= kl_fctr*self.variational_prior.KL_divergence(X) + + current_values['meangrad'], current_values['vargrad'] = self.kern.gradients_qX_expectations( + variational_posterior=X, + Z=Z, dL_dpsi0=grad_dict['dL_dpsi0'], + dL_dpsi1=grad_dict['dL_dpsi1'], + dL_dpsi2=grad_dict['dL_dpsi2']) + + # Subsetting Variational Posterior objects, makes the gradients + # empty. We need them to be 0 though: + X.mean.gradient[:] = 0 + X.variance.gradient[:] = 0 + + self.variational_prior.update_gradients_KL(X) + if self.missing_data: + current_values['meangrad'] += kl_fctr*X.mean.gradient/d + current_values['vargrad'] += kl_fctr*X.variance.gradient/d + else: + current_values['meangrad'] += kl_fctr*X.mean.gradient + current_values['vargrad'] += kl_fctr*X.variance.gradient + + if subset_indices is not None: + value_indices['meangrad'] = subset_indices['samples'] + value_indices['vargrad'] = subset_indices['samples'] + return posterior, log_marginal_likelihood, grad_dict, current_values, value_indices + + def _outer_values_update(self, full_values): + """ + Here you put the values, which were collected before in the right places. + E.g. set the gradients of parameters, etc. + """ + super(BayesianGPLVMMiniBatch, self)._outer_values_update(full_values) + self.X.mean.gradient = full_values['meangrad'] + self.X.variance.gradient = full_values['vargrad'] + + def _outer_init_full_values(self): + return dict(meangrad=np.zeros(self.X.mean.shape), + vargrad=np.zeros(self.X.variance.shape)) + + def parameters_changed(self): + super(BayesianGPLVMMiniBatch,self).parameters_changed() + if isinstance(self.inference_method, VarDTC_minibatch): + return + + #super(BayesianGPLVM, self).parameters_changed() + #self._log_marginal_likelihood -= self.variational_prior.KL_divergence(self.X) + + #self.X.mean.gradient, self.X.variance.gradient = self.kern.gradients_qX_expectations(variational_posterior=self.X, Z=self.Z, dL_dpsi0=self.grad_dict['dL_dpsi0'], dL_dpsi1=self.grad_dict['dL_dpsi1'], dL_dpsi2=self.grad_dict['dL_dpsi2']) + + # This is testing code ------------------------- +# i = np.random.randint(self.X.shape[0]) +# X_ = self.X.mean +# which = np.sqrt(((X_ - X_[i:i+1])**2).sum(1)).argsort()>(max(0, self.X.shape[0]-51)) +# _, _, grad_dict = self.inference_method.inference(self.kern, self.X[which], self.Z, self.likelihood, self.Y[which], self.Y_metadata) +# grad = self.kern.gradients_qX_expectations(variational_posterior=self.X[which], Z=self.Z, dL_dpsi0=grad_dict['dL_dpsi0'], dL_dpsi1=grad_dict['dL_dpsi1'], dL_dpsi2=grad_dict['dL_dpsi2']) +# +# self.X.mean.gradient[:] = 0 +# self.X.variance.gradient[:] = 0 +# self.X.mean.gradient[which] = grad[0] +# self.X.variance.gradient[which] = grad[1] + + # update for the KL divergence +# self.variational_prior.update_gradients_KL(self.X, which) + # ----------------------------------------------- + + # update for the KL divergence + #self.variational_prior.update_gradients_KL(self.X) + + def plot_latent(self, labels=None, which_indices=None, + resolution=50, ax=None, marker='o', s=40, + fignum=None, plot_inducing=True, legend=True, + plot_limits=None, + aspect='auto', updates=False, predict_kwargs={}, imshow_kwargs={}): + import sys + assert "matplotlib" in sys.modules, "matplotlib package has not been imported." + from ..plotting.matplot_dep import dim_reduction_plots + + return dim_reduction_plots.plot_latent(self, labels, which_indices, + resolution, ax, marker, s, + fignum, plot_inducing, legend, + plot_limits, aspect, updates, predict_kwargs, imshow_kwargs) + + def do_test_latents(self, Y): + """ + Compute the latent representation for a set of new points Y + + Notes: + This will only work with a univariate Gaussian likelihood (for now) + """ + N_test = Y.shape[0] + input_dim = self.Z.shape[1] + + means = np.zeros((N_test, input_dim)) + covars = np.zeros((N_test, input_dim)) + + dpsi0 = -0.5 * self.input_dim / self.likelihood.variance + dpsi2 = self.grad_dict['dL_dpsi2'][0][None, :, :] # TODO: this may change if we ignore het. likelihoods + V = Y/self.likelihood.variance + + #compute CPsi1V + #if self.Cpsi1V is None: + # psi1V = np.dot(self.psi1.T, self.likelihood.V) + # tmp, _ = linalg.dtrtrs(self._Lm, np.asfortranarray(psi1V), lower=1, trans=0) + # tmp, _ = linalg.dpotrs(self.LB, tmp, lower=1) + # self.Cpsi1V, _ = linalg.dtrtrs(self._Lm, tmp, lower=1, trans=1) + + dpsi1 = np.dot(self.posterior.woodbury_vector, V.T) + + #start = np.zeros(self.input_dim * 2) + + + from scipy.optimize import minimize + + for n, dpsi1_n in enumerate(dpsi1.T[:, :, None]): + args = (input_dim, self.kern.copy(), self.Z, dpsi0, dpsi1_n.T, dpsi2) + res = minimize(latent_cost_and_grad, jac=True, x0=np.hstack((means[n], covars[n])), args=args, method='BFGS') + xopt = res.x + mu, log_S = xopt.reshape(2, 1, -1) + means[n] = mu[0].copy() + covars[n] = np.exp(log_S[0]).copy() + + X = NormalPosterior(means, covars) + + return X + + def dmu_dX(self, Xnew): + """ + Calculate the gradient of the prediction at Xnew w.r.t Xnew. + """ + dmu_dX = np.zeros_like(Xnew) + for i in range(self.Z.shape[0]): + dmu_dX += self.kern.gradients_X(self.grad_dict['dL_dpsi1'][i:i + 1, :], Xnew, self.Z[i:i + 1, :]) + return dmu_dX + + def dmu_dXnew(self, Xnew): + """ + Individual gradient of prediction at Xnew w.r.t. each sample in Xnew + """ + gradients_X = np.zeros((Xnew.shape[0], self.num_inducing)) + ones = np.ones((1, 1)) + for i in range(self.Z.shape[0]): + gradients_X[:, i] = self.kern.gradients_X(ones, Xnew, self.Z[i:i + 1, :]).sum(-1) + return np.dot(gradients_X, self.grad_dict['dL_dpsi1']) + + def plot_steepest_gradient_map(self, *args, ** kwargs): + """ + See GPy.plotting.matplot_dep.dim_reduction_plots.plot_steepest_gradient_map + """ + import sys + assert "matplotlib" in sys.modules, "matplotlib package has not been imported." + from ..plotting.matplot_dep import dim_reduction_plots + + return dim_reduction_plots.plot_steepest_gradient_map(self,*args,**kwargs) + + +def latent_cost_and_grad(mu_S, input_dim, kern, Z, dL_dpsi0, dL_dpsi1, dL_dpsi2): + """ + objective function for fitting the latent variables for test points + (negative log-likelihood: should be minimised!) + """ + mu = mu_S[:input_dim][None] + log_S = mu_S[input_dim:][None] + S = np.exp(log_S) + + X = NormalPosterior(mu, S) + + psi0 = kern.psi0(Z, X) + psi1 = kern.psi1(Z, X) + psi2 = kern.psi2(Z, X) + + lik = dL_dpsi0 * psi0.sum() + np.einsum('ij,kj->...', dL_dpsi1, psi1) + np.einsum('ijk,lkj->...', dL_dpsi2, psi2) - 0.5 * np.sum(np.square(mu) + S) + 0.5 * np.sum(log_S) + + dLdmu, dLdS = kern.gradients_qX_expectations(dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, X) + dmu = dLdmu - mu + # dS = S0 + S1 + S2 -0.5 + .5/S + dlnS = S * (dLdS - 0.5) + .5 + + return -lik, -np.hstack((dmu.flatten(), dlnS.flatten())) diff --git a/GPy/models/sparse_gp_minibatch.py b/GPy/models/sparse_gp_minibatch.py new file mode 100644 index 00000000..00929462 --- /dev/null +++ b/GPy/models/sparse_gp_minibatch.py @@ -0,0 +1,347 @@ +# Copyright (c) 2012, GPy authors (see AUTHORS.txt). +# Licensed under the BSD 3-clause license (see LICENSE.txt) + +import numpy as np +from ..core.parameterization.param import Param +from ..core.gp import GP +from ..inference.latent_function_inference import var_dtc +from .. import likelihoods +from ..core.parameterization.variational import VariationalPosterior + +import logging +from GPy.inference.latent_function_inference.posterior import Posterior +from GPy.inference.optimization.stochastics import SparseGPStochastics,\ + SparseGPMissing +#no stochastics.py file added! from GPy.inference.optimization.stochastics import SparseGPStochastics,\ + #SparseGPMissing +logger = logging.getLogger("sparse gp") + +class SparseGPMiniBatch(GP): + """ + A general purpose Sparse GP model +''' +Created on 3 Nov 2014 + +@author: maxz +''' + + This model allows (approximate) inference using variational DTC or FITC + (Gaussian likelihoods) as well as non-conjugate sparse methods based on + these. + + :param X: inputs + :type X: np.ndarray (num_data x input_dim) + :param likelihood: a likelihood instance, containing the observed data + :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 (num_data x input_dim) | None + :param Z: inducing inputs + :type Z: np.ndarray (num_inducing x input_dim) + :param num_inducing: Number of inducing points (optional, default 10. Ignored if Z is not None) + :type num_inducing: int + + """ + + def __init__(self, X, Y, Z, kernel, likelihood, inference_method=None, + name='sparse gp', Y_metadata=None, normalizer=False, + missing_data=False, stochastic=False, batchsize=1): + #pick a sensible inference method + if inference_method is None: + if isinstance(likelihood, likelihoods.Gaussian): + inference_method = var_dtc.VarDTC(limit=1 if not self.missing_data else Y.shape[1]) + else: + #inference_method = ?? + raise NotImplementedError, "what to do what to do?" + print "defaulting to ", inference_method, "for latent function inference" + + self.Z = Param('inducing inputs', Z) + self.num_inducing = Z.shape[0] + + GP.__init__(self, X, Y, kernel, likelihood, inference_method=inference_method, name=name, Y_metadata=Y_metadata, normalizer=normalizer) + self.missing_data = missing_data + + if stochastic and missing_data: + self.missing_data = True + self.ninan = ~np.isnan(Y) + self.stochastics = SparseGPStochastics(self, batchsize) + elif stochastic and not missing_data: + self.missing_data = False + self.stochastics = SparseGPStochastics(self, batchsize) + elif missing_data: + self.missing_data = True + self.ninan = ~np.isnan(Y) + self.stochastics = SparseGPMissing(self) + else: + self.stochastics = False + + logger.info("Adding Z as parameter") + self.link_parameter(self.Z, index=0) + if self.missing_data: + self.Ylist = [] + overall = self.Y_normalized.shape[1] + m_f = lambda i: "Precomputing Y for missing data: {: >7.2%}".format(float(i+1)/overall) + message = m_f(-1) + print message, + for d in xrange(overall): + self.Ylist.append(self.Y_normalized[self.ninan[:, d], d][:, None]) + print ' '*(len(message)+1) + '\r', + message = m_f(d) + print message, + print '' + + self.posterior = None + + def has_uncertain_inputs(self): + return isinstance(self.X, VariationalPosterior) + + def _inner_parameters_changed(self, kern, X, Z, likelihood, Y, Y_metadata, Lm=None, dL_dKmm=None, subset_indices=None): + """ + This is the standard part, which usually belongs in parameters_changed. + + For automatic handling of subsampling (such as missing_data, stochastics etc.), we need to put this into an inner + loop, in order to ensure a different handling of gradients etc of different + subsets of data. + + The dict in current_values will be passed aroung as current_values for + the rest of the algorithm, so this is the place to store current values, + such as subsets etc, if necessary. + + If Lm and dL_dKmm can be precomputed (or only need to be computed once) + pass them in here, so they will be passed to the inference_method. + + subset_indices is a dictionary of indices. you can put the indices however you + like them into this dictionary for inner use of the indices inside the + algorithm. + """ + try: + posterior, log_marginal_likelihood, grad_dict = self.inference_method.inference(kern, X, Z, likelihood, Y, Y_metadata, Lm=Lm, dL_dKmm=None) + except: + posterior, log_marginal_likelihood, grad_dict = self.inference_method.inference(kern, X, Z, likelihood, Y, Y_metadata) + current_values = {} + likelihood.update_gradients(grad_dict['dL_dthetaL']) + current_values['likgrad'] = likelihood.gradient.copy() + if subset_indices is None: + subset_indices = {} + if isinstance(X, VariationalPosterior): + #gradients wrt kernel + dL_dKmm = grad_dict['dL_dKmm'] + kern.update_gradients_full(dL_dKmm, Z, None) + current_values['kerngrad'] = kern.gradient.copy() + kern.update_gradients_expectations(variational_posterior=X, + Z=Z, + dL_dpsi0=grad_dict['dL_dpsi0'], + dL_dpsi1=grad_dict['dL_dpsi1'], + dL_dpsi2=grad_dict['dL_dpsi2']) + current_values['kerngrad'] += kern.gradient + + #gradients wrt Z + current_values['Zgrad'] = kern.gradients_X(dL_dKmm, Z) + current_values['Zgrad'] += kern.gradients_Z_expectations( + grad_dict['dL_dpsi0'], + grad_dict['dL_dpsi1'], + grad_dict['dL_dpsi2'], + Z=Z, + variational_posterior=X) + else: + #gradients wrt kernel + kern.update_gradients_diag(grad_dict['dL_dKdiag'], X) + current_values['kerngrad'] = kern.gradient.copy() + kern.update_gradients_full(grad_dict['dL_dKnm'], X, Z) + current_values['kerngrad'] += kern.gradient + kern.update_gradients_full(grad_dict['dL_dKmm'], Z, None) + current_values['kerngrad'] += kern.gradient + #gradients wrt Z + current_values['Zgrad'] = kern.gradients_X(grad_dict['dL_dKmm'], Z) + current_values['Zgrad'] += kern.gradients_X(grad_dict['dL_dKnm'].T, Z, X) + return posterior, log_marginal_likelihood, grad_dict, current_values, subset_indices + + def _inner_take_over_or_update(self, full_values=None, current_values=None, value_indices=None): + """ + This is for automatic updates of values in the inner loop of missing + data handling. Both arguments are dictionaries and the values in + full_values will be updated by the current_gradients. + + If a key from current_values does not exist in full_values, it will be + initialized to the value in current_values. + + If there is indices needed for the update, value_indices can be used for + that. If value_indices has the same key, as current_values, the update + in full_values will be indexed by the indices in value_indices. + + grads: + dictionary of standing gradients (you will have to carefully make sure, that + the ordering is right!). The values in here will be updated such that + full_values[key] += current_values[key] forall key in full_gradients.keys() + + gradients: + dictionary of gradients in the current set of parameters. + + value_indices: + dictionary holding indices for the update in full_values. + if the key exists the update rule is:def df(x): + full_values[key][value_indices[key]] += current_values[key] + """ + for key in current_values.keys(): + if value_indices is not None and value_indices.has_key(key): + index = value_indices[key] + else: + index = slice(None) + if full_values.has_key(key): + full_values[key][index] += current_values[key] + else: + full_values[key] = current_values[key] + + def _inner_values_update(self, current_values): + """ + This exists if there is more to do with the current values. + It will be called allways in the inner loop, so that + you can do additional inner updates for the inside of the missing data + loop etc. This can also be used for stochastic updates, when only working on + one dimension of the output. + """ + pass + + def _outer_values_update(self, full_values): + """ + Here you put the values, which were collected before in the right places. + E.g. set the gradients of parameters, etc. + """ + self.likelihood.gradient = full_values['likgrad'] + self.kern.gradient = full_values['kerngrad'] + self.Z.gradient = full_values['Zgrad'] + + def _outer_init_full_values(self): + """ + If full_values has indices in values_indices, we might want to initialize + the full_values differently, so that subsetting is possible. + + Here you can initialize the full_values for the values needed. + + Keep in mind, that if a key does not exist in full_values when updating + values, it will be set (so e.g. for Z there is no need to initialize Zgrad, + as there is no subsetting needed. For X in BGPLVM on the other hand we probably need + to initialize the gradients for the mean and the variance in order to + have the full gradient for indexing) + """ + return {} + + def _outer_loop_for_missing_data(self): + Lm = None + dL_dKmm = None + + self._log_marginal_likelihood = 0 + self.full_values = self._outer_init_full_values() + + if self.posterior is None: + woodbury_inv = np.zeros((self.num_inducing, self.num_inducing, self.output_dim)) + woodbury_vector = np.zeros((self.num_inducing, self.output_dim)) + else: + woodbury_inv = self.posterior._woodbury_inv + woodbury_vector = self.posterior._woodbury_vector + + if not self.stochastics: + m_f = lambda i: "Inference with missing_data: {: >7.2%}".format(float(i+1)/self.output_dim) + message = m_f(-1) + print message, + + for d in self.stochastics.d: + ninan = self.ninan[:, d] + + if not self.stochastics: + print ' '*(len(message)) + '\r', + message = m_f(d) + print message, + + posterior, log_marginal_likelihood, \ + grad_dict, current_values, value_indices = self._inner_parameters_changed( + self.kern, self.X[ninan], + self.Z, self.likelihood, + self.Ylist[d], self.Y_metadata, + Lm, dL_dKmm, + subset_indices=dict(outputs=d, samples=ninan)) + + self._inner_take_over_or_update(self.full_values, current_values, value_indices) + self._inner_values_update(current_values) + + Lm = posterior.K_chol + dL_dKmm = grad_dict['dL_dKmm'] + woodbury_inv[:, :, d] = posterior.woodbury_inv + woodbury_vector[:, d:d+1] = posterior.woodbury_vector + self._log_marginal_likelihood += log_marginal_likelihood + if not self.stochastics: + print '' + + if self.posterior is None: + self.posterior = Posterior(woodbury_inv=woodbury_inv, woodbury_vector=woodbury_vector, + K=posterior._K, mean=None, cov=None, K_chol=posterior.K_chol) + self._outer_values_update(self.full_values) + + def _outer_loop_without_missing_data(self): + self._log_marginal_likelihood = 0 + + if self.posterior is None: + woodbury_inv = np.zeros((self.num_inducing, self.num_inducing, self.output_dim)) + woodbury_vector = np.zeros((self.num_inducing, self.output_dim)) + else: + woodbury_inv = self.posterior._woodbury_inv + woodbury_vector = self.posterior._woodbury_vector + + d = self.stochastics.d + posterior, log_marginal_likelihood, \ + grad_dict, self.full_values, _ = self._inner_parameters_changed( + self.kern, self.X, + self.Z, self.likelihood, + self.Y_normalized[:, d], self.Y_metadata) + self.grad_dict = grad_dict + + self._log_marginal_likelihood += log_marginal_likelihood + + self._outer_values_update(self.full_values) + + woodbury_inv[:, :, d] = posterior.woodbury_inv[:, :, None] + woodbury_vector[:, d] = posterior.woodbury_vector + if self.posterior is None: + self.posterior = Posterior(woodbury_inv=woodbury_inv, woodbury_vector=woodbury_vector, + K=posterior._K, mean=None, cov=None, K_chol=posterior.K_chol) + + def parameters_changed(self): + if self.missing_data: + self._outer_loop_for_missing_data() + elif self.stochastics: + self._outer_loop_without_missing_data() + else: + self.posterior, self._log_marginal_likelihood, self.grad_dict, self.full_values, _ = self._inner_parameters_changed(self.kern, self.X, self.Z, self.likelihood, self.Y_normalized, self.Y_metadata) + self._outer_values_update(self.full_values) + + def _raw_predict(self, Xnew, full_cov=False, kern=None): + """ + Make a prediction for the latent function values + """ + + if kern is None: kern = self.kern + + if not isinstance(Xnew, VariationalPosterior): + Kx = kern.K(self.Z, Xnew) + mu = np.dot(Kx.T, self.posterior.woodbury_vector) + if full_cov: + Kxx = kern.K(Xnew) + if self.posterior.woodbury_inv.ndim == 2: + var = Kxx - np.dot(Kx.T, np.dot(self.posterior.woodbury_inv, Kx)) + elif self.posterior.woodbury_inv.ndim == 3: + var = Kxx[:,:,None] - np.tensordot(np.dot(np.atleast_3d(self.posterior.woodbury_inv).T, Kx).T, Kx, [1,0]).swapaxes(1,2) + var = var + else: + Kxx = kern.Kdiag(Xnew) + var = (Kxx - np.sum(np.dot(np.atleast_3d(self.posterior.woodbury_inv).T, Kx) * Kx[None,:,:], 1)).T + else: + Kx = kern.psi1(self.Z, Xnew) + mu = np.dot(Kx, self.posterior.woodbury_vector) + if full_cov: + raise NotImplementedError, "TODO" + else: + Kxx = kern.psi0(self.Z, Xnew) + psi2 = kern.psi2(self.Z, Xnew) + var = Kxx - np.sum(np.sum(psi2 * Kmmi_LmiBLmi[None, :, :], 1), 1) + return mu, var diff --git a/GPy/testing/model_tests.py b/GPy/testing/model_tests.py index d05e33de..55dfafec 100644 --- a/GPy/testing/model_tests.py +++ b/GPy/testing/model_tests.py @@ -112,7 +112,7 @@ class MiscTests(unittest.TestCase): def test_missing_data(self): from GPy import kern - from GPy.models import BayesianGPLVM + from GPy.models.bayesian_gplvm_minibatch import BayesianGPLVMMiniBatch from GPy.examples.dimensionality_reduction import _simulate_matern D1, D2, D3, N, num_inducing, Q = 13, 5, 8, 400, 3, 4 @@ -124,12 +124,12 @@ class MiscTests(unittest.TestCase): Ymissing[inan] = np.nan k = kern.Linear(Q, ARD=True) + kern.White(Q, np.exp(-2)) # + kern.bias(Q) - m = BayesianGPLVM(Ymissing, Q, init="random", num_inducing=num_inducing, + m = BayesianGPLVMMiniBatch(Ymissing, Q, init="random", num_inducing=num_inducing, kernel=k, missing_data=True) assert(m.checkgrad()) k = kern.RBF(Q, ARD=True) + kern.White(Q, np.exp(-2)) # + kern.bias(Q) - m = BayesianGPLVM(Ymissing, Q, init="random", num_inducing=num_inducing, + m = BayesianGPLVMMiniBatch(Ymissing, Q, init="random", num_inducing=num_inducing, kernel=k, missing_data=True) assert(m.checkgrad()) From 02c903c4eb732399bf5aaa98abc6e6a2fdd8db8b Mon Sep 17 00:00:00 2001 From: Alan Saul Date: Mon, 3 Nov 2014 13:36:26 +0000 Subject: [PATCH 2/5] Tidied up laplace warnings --- .../latent_function_inference/laplace.py | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/GPy/inference/latent_function_inference/laplace.py b/GPy/inference/latent_function_inference/laplace.py index 01d2f997..2c741b9d 100644 --- a/GPy/inference/latent_function_inference/laplace.py +++ b/GPy/inference/latent_function_inference/laplace.py @@ -14,6 +14,9 @@ import numpy as np from ...util.linalg import mdot, jitchol, dpotrs, dtrtrs, dpotri, symmetrify, pdinv from posterior import Posterior import warnings +def warning_on_one_line(message, category, filename, lineno, file=None, line=None): + return ' %s:%s: %s:%s\n' % (filename, lineno, category.__name__, message) +warnings.formatwarning = warning_on_one_line from scipy import optimize from . import LatentFunctionInference @@ -29,8 +32,11 @@ class Laplace(LatentFunctionInference): """ self._mode_finding_tolerance = 1e-7 - self._mode_finding_max_iter = 40 - self.bad_fhat = True + self._mode_finding_max_iter = 60 + self.bad_fhat = False + #Store whether it is the first run of the inference so that we can choose whether we need + #to calculate things or reuse old variables + self.first_run = True self._previous_Ki_fhat = None def inference(self, kern, X, likelihood, Y, Y_metadata=None): @@ -42,8 +48,9 @@ class Laplace(LatentFunctionInference): K = kern.K(X) #Find mode - if self.bad_fhat: + if self.bad_fhat or self.first_run: Ki_f_init = np.zeros_like(Y) + first_run = False else: Ki_f_init = self._previous_Ki_fhat @@ -123,11 +130,11 @@ class Laplace(LatentFunctionInference): #Warn of bad fits if difference > self._mode_finding_tolerance: if not self.bad_fhat: - warnings.warn("Not perfect f_hat fit difference: {}".format(difference)) + warnings.warn("Not perfect mode found (f_hat). difference: {}, iteration: {} out of max {}".format(difference, iteration, self._mode_finding_max_iter)) self.bad_fhat = True elif self.bad_fhat: self.bad_fhat = False - warnings.warn("f_hat now fine again") + warnings.warn("f_hat now fine again. difference: {}, iteration: {} out of max {}".format(difference, iteration, self._mode_finding_max_iter)) return f, Ki_f From 1385c3b2711277034e654874090da277cc25a277 Mon Sep 17 00:00:00 2001 From: Alan Saul Date: Mon, 3 Nov 2014 13:36:46 +0000 Subject: [PATCH 3/5] Removed pod dependency for pickle tests --- GPy/examples/regression.py | 2 +- GPy/testing/pickle_tests.py | 17 +++++++++++------ 2 files changed, 12 insertions(+), 7 deletions(-) diff --git a/GPy/examples/regression.py b/GPy/examples/regression.py index 8e481b02..6923b20d 100644 --- a/GPy/examples/regression.py +++ b/GPy/examples/regression.py @@ -468,7 +468,7 @@ def uncertain_inputs_sparse_regression(max_iters=200, optimize=True, plot=True): k = GPy.kern.RBF(1) # create simple GP Model - no input uncertainty on this one - m = GPy.models.SparseGPRegression(X, Y, kernel=GPy.kern.RBF(1), Z=Z) + m = GPy.models.SparseGPRegression(X, Y, kernel=k, Z=Z) if optimize: m.optimize('scg', messages=1, max_iters=max_iters) diff --git a/GPy/testing/pickle_tests.py b/GPy/testing/pickle_tests.py index dfabe54e..b541d21e 100644 --- a/GPy/testing/pickle_tests.py +++ b/GPy/testing/pickle_tests.py @@ -17,10 +17,15 @@ from GPy.kern._src.rbf import RBF from GPy.kern._src.linear import Linear from GPy.kern._src.static import Bias, White from GPy.examples.dimensionality_reduction import mrd_simulation -from GPy.examples.regression import toy_rbf_1d_50 from GPy.core.parameterization.variational import NormalPosterior from GPy.models.gp_regression import GPRegression +def toy_model(): + X = np.linspace(0,1,50)[:, None] + Y = np.sin(X) + m = GPRegression(X=X, Y=Y) + return m + class ListDictTestCase(unittest.TestCase): def assertListDictEquals(self, d1, d2, msg=None): for k,v in d1.iteritems(): @@ -105,7 +110,7 @@ class Test(ListDictTestCase): self.assertSequenceEqual(str(par), str(pcopy)) def test_model(self): - par = toy_rbf_1d_50(optimize=0, plot=0) + par = toy_model() pcopy = par.copy() self.assertListEqual(par.param_array.tolist(), pcopy.param_array.tolist()) np.testing.assert_allclose(par.gradient_full, pcopy.gradient_full) @@ -124,7 +129,7 @@ class Test(ListDictTestCase): self.assert_(pcopy.checkgrad()) def test_modelrecreation(self): - par = toy_rbf_1d_50(optimize=0, plot=0) + par = toy_model() pcopy = GPRegression(par.X.copy(), par.Y.copy(), kernel=par.kern.copy()) np.testing.assert_allclose(par.param_array, pcopy.param_array) np.testing.assert_allclose(par.gradient_full, pcopy.gradient_full) @@ -135,7 +140,7 @@ class Test(ListDictTestCase): self.assert_(np.any(pcopy.gradient!=0.0)) pcopy.optimize('bfgs') par.optimize('bfgs') - np.testing.assert_allclose(pcopy.param_array, par.param_array, atol=.001) + np.testing.assert_allclose(pcopy.param_array, par.param_array, atol=1e-6) with tempfile.TemporaryFile('w+b') as f: par.pickle(f) f.seek(0) @@ -193,7 +198,7 @@ class Test(ListDictTestCase): @unittest.skip def test_add_observer(self): - par = toy_rbf_1d_50(optimize=0, plot=0) + par = toy_model() par.name = "original" par.count = 0 par.add_observer(self, self._callback, 1) @@ -211,4 +216,4 @@ class Test(ListDictTestCase): if __name__ == "__main__": #import sys;sys.argv = ['', 'Test.test_parameter_index_operations'] - unittest.main() \ No newline at end of file + unittest.main() From d0a5420f2f3e27931065c4f79c0146135d4620bf Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Mon, 3 Nov 2014 14:19:03 +0000 Subject: [PATCH 4/5] [minbatch var dtc] adjustments to bgplvm minibatch --- GPy/models/bayesian_gplvm.py | 64 ++++++-------------------- GPy/models/bayesian_gplvm_minibatch.py | 4 +- 2 files changed, 17 insertions(+), 51 deletions(-) diff --git a/GPy/models/bayesian_gplvm.py b/GPy/models/bayesian_gplvm.py index c4d0b6e9..73fd6f8f 100644 --- a/GPy/models/bayesian_gplvm.py +++ b/GPy/models/bayesian_gplvm.py @@ -75,8 +75,7 @@ class BayesianGPLVM(SparseGP_MPI): name=name, inference_method=inference_method, normalizer=normalizer, mpi_comm=mpi_comm, variational_prior=self.variational_prior, - missing_data=missing_data, stochastic=stochastic, - batchsize=batchsize) + ) def set_X_gradients(self, X, X_grad): """Set the gradients of the posterior distribution of X in its specific form.""" @@ -86,55 +85,22 @@ class BayesianGPLVM(SparseGP_MPI): """Get the gradients of the posterior distribution of X in its specific form.""" return X.mean.gradient, X.variance.gradient - def _inner_parameters_changed(self, kern, X, Z, likelihood, Y, Y_metadata, Lm=None, dL_dKmm=None, subset_indices=None): - posterior, log_marginal_likelihood, grad_dict, current_values, value_indices = super(BayesianGPLVM, self)._inner_parameters_changed(kern, X, Z, likelihood, Y, Y_metadata, Lm=Lm, dL_dKmm=dL_dKmm, subset_indices=subset_indices) - - kl_fctr = 1. - if self.missing_data: - d = self.output_dim - log_marginal_likelihood -= kl_fctr*self.variational_prior.KL_divergence(X)/d - else: - log_marginal_likelihood -= kl_fctr*self.variational_prior.KL_divergence(X) - - current_values['meangrad'], current_values['vargrad'] = self.kern.gradients_qX_expectations( - variational_posterior=X, - Z=Z, dL_dpsi0=grad_dict['dL_dpsi0'], - dL_dpsi1=grad_dict['dL_dpsi1'], - dL_dpsi2=grad_dict['dL_dpsi2']) - - # Subsetting Variational Posterior objects, makes the gradients - # empty. We need them to be 0 though: - X.mean.gradient[:] = 0 - X.variance.gradient[:] = 0 - - self.variational_prior.update_gradients_KL(X) - if self.missing_data: - current_values['meangrad'] += kl_fctr*X.mean.gradient/d - current_values['vargrad'] += kl_fctr*X.variance.gradient/d - else: - current_values['meangrad'] += kl_fctr*X.mean.gradient - current_values['vargrad'] += kl_fctr*X.variance.gradient - - if subset_indices is not None: - value_indices['meangrad'] = subset_indices['samples'] - value_indices['vargrad'] = subset_indices['samples'] - return posterior, log_marginal_likelihood, grad_dict, current_values, value_indices - - def _outer_values_update(self, full_values): - """ - Here you put the values, which were collected before in the right places. - E.g. set the gradients of parameters, etc. - """ - super(BayesianGPLVM, self)._outer_values_update(full_values) - self.X.mean.gradient = full_values['meangrad'] - self.X.variance.gradient = full_values['vargrad'] - - def _outer_init_full_values(self): - return dict(meangrad=np.zeros(self.X.mean.shape), - vargrad=np.zeros(self.X.variance.shape)) - def parameters_changed(self): super(BayesianGPLVM,self).parameters_changed() + + kl_fctr = 1. + self._log_marginal_likelihood -= kl_fctr*self.variational_prior.KL_divergence(self.X) + + self.X.mean.gradient, self.X.variance.gradient = self.kern.gradients_qX_expectations( + variational_posterior=self.X, + Z=self.Z, + dL_dpsi0=self.grad_dict['dL_dpsi0'], + dL_dpsi1=self.grad_dict['dL_dpsi1'], + dL_dpsi2=self.grad_dict['dL_dpsi2']) + + self.variational_prior.update_gradients_KL(self.X) + + if isinstance(self.inference_method, VarDTC_minibatch): return diff --git a/GPy/models/bayesian_gplvm_minibatch.py b/GPy/models/bayesian_gplvm_minibatch.py index 56b46244..7100fa72 100644 --- a/GPy/models/bayesian_gplvm_minibatch.py +++ b/GPy/models/bayesian_gplvm_minibatch.py @@ -26,8 +26,6 @@ class BayesianGPLVMMiniBatch(SparseGPMiniBatch): Z=None, kernel=None, inference_method=None, likelihood=None, name='bayesian gplvm', normalizer=None, missing_data=False, stochastic=False, batchsize=1): - self.__IN_OPTIMIZATION__ = False - self.logger = logging.getLogger(self.__class__.__name__) if X is None: from ..util.initialization import initialize_latent @@ -70,6 +68,8 @@ class BayesianGPLVMMiniBatch(SparseGPMiniBatch): normalizer=normalizer, missing_data=missing_data, stochastic=stochastic, batchsize=batchsize) + self.X = X + self.link_parameter(self.X, 0) def set_X_gradients(self, X, X_grad): """Set the gradients of the posterior distribution of X in its specific form.""" From 6242a75f6a28b41fbe699a1210442b5d7354c12f Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Mon, 3 Nov 2014 14:19:38 +0000 Subject: [PATCH 5/5] [parameterized] handle updates inside init --- GPy/core/parameterization/index_operations.py | 20 +++++++++---------- GPy/core/parameterization/parameter_core.py | 6 +++--- GPy/core/parameterization/parameterized.py | 15 ++++++++++---- 3 files changed, 23 insertions(+), 18 deletions(-) diff --git a/GPy/core/parameterization/index_operations.py b/GPy/core/parameterization/index_operations.py index 78c0d2b9..ddd689ed 100644 --- a/GPy/core/parameterization/index_operations.py +++ b/GPy/core/parameterization/index_operations.py @@ -26,39 +26,39 @@ class ParameterIndexOperations(object): This object wraps a dictionary, whos keys are _operations_ that we'd like to apply to a parameter array, and whose values are np integer arrays which index the parameter array appropriately. - + A model instance will contain one instance of this class for each thing that needs indexing (i.e. constraints, ties and priors). Parameters within the model constain instances of the ParameterIndexOperationsView class, which can map from a 'local' index (starting 0) to this global index. - + Here's an illustration: - + #======================================================================= model : 0 1 2 3 4 5 6 7 8 9 key1: 4 5 key2: 7 8 - + param1: 0 1 2 3 4 5 key1: 2 3 key2: 5 - + param2: 0 1 2 3 4 key1: 0 key2: 2 3 #======================================================================= - + The views of this global index have a subset of the keys in this global (model) index. - + Adding a new key (e.g. a constraint) to a view will cause the view to pass the new key to the global index, along with the local index and an offset. This global index then stores the key and the appropriate global index (which can be seen by the view). - + See also: ParameterIndexOperationsView - + """ _offset = 0 def __init__(self, constraints=None): @@ -221,8 +221,6 @@ class ParameterIndexOperationsView(object): def shift_left(self, start, size): self._param_index_ops.shift_left(start+self._offset, size) - self._offset -= size - self._size -= size def clear(self): for i, ind in self.items(): diff --git a/GPy/core/parameterization/parameter_core.py b/GPy/core/parameterization/parameter_core.py index b1e0d8d2..115acf3a 100644 --- a/GPy/core/parameterization/parameter_core.py +++ b/GPy/core/parameterization/parameter_core.py @@ -18,7 +18,7 @@ import numpy as np import re import logging -__updated__ = '2014-10-28' +__updated__ = '2014-11-03' class HierarchyError(Exception): """ @@ -924,7 +924,7 @@ class Parameterizable(OptimizationHandlable): !WARNING!: setting the parameter array MUST always be done in memory: m.param_array[:] = m_copy.param_array """ - if self.__dict__.get('_param_array_', None) is None: + if (self.__dict__.get('_param_array_', None) is None) or (self._param_array_.size != self.size): self._param_array_ = np.empty(self.size, dtype=np.float64) return self._param_array_ @@ -1002,7 +1002,7 @@ class Parameterizable(OptimizationHandlable): #========================================================================= @property def gradient(self): - if self.__dict__.get('_gradient_array_', None) is None: + if (self.__dict__.get('_gradient_array_', None) is None) or self._gradient_array_.size != self.size: self._gradient_array_ = np.empty(self.size, dtype=np.float64) return self._gradient_array_ diff --git a/GPy/core/parameterization/parameterized.py b/GPy/core/parameterization/parameterized.py index 0bd72be4..2c91b235 100644 --- a/GPy/core/parameterization/parameterized.py +++ b/GPy/core/parameterization/parameterized.py @@ -9,6 +9,7 @@ from param import ParamConcatenation from parameter_core import HierarchyError, Parameterizable, adjust_name_for_printing import logging +from GPy.core.parameterization.index_operations import ParameterIndexOperationsView logger = logging.getLogger("parameters changed meta") class ParametersChangedMeta(type): @@ -20,7 +21,7 @@ class ParametersChangedMeta(type): self._in_init_ = False logger.debug("connecting parameters") self._highest_parent_._connect_parameters() - self._highest_parent_._notify_parent_change() + #self._highest_parent_._notify_parent_change() self._highest_parent_._connect_fixes() logger.debug("calling parameters changed") self.parameters_changed() @@ -140,6 +141,8 @@ class Parameterized(Parameterizable): self.priors.shift_right(start, param.size) self.constraints.update(param.constraints, self.size) self.priors.update(param.priors, self.size) + param._parent_ = self + param._parent_index_ = len(self.parameters) self.parameters.append(param) else: start = sum(p.size for p in self.parameters[:index]) @@ -147,19 +150,23 @@ class Parameterized(Parameterizable): self.priors.shift_right(start, param.size) self.constraints.update(param.constraints, start) self.priors.update(param.priors, start) + param._parent_ = self + param._parent_index_ = index if index>=0 else len(self.parameters[:index]) + for p in self.parameters[index:]: + p._parent_index_ += 1 self.parameters.insert(index, param) - self._notify_parent_change() param.add_observer(self, self._pass_through_notify_observers, -np.inf) parent = self while parent is not None: parent.size += param.size parent = parent._parent_ + self._notify_parent_change() if not self._in_init_: - self._connect_parameters() - self._notify_parent_change() + #self._connect_parameters() + #self._notify_parent_change() self._highest_parent_._connect_parameters(ignore_added_names=_ignore_added_names) self._highest_parent_._notify_parent_change()