From b20beaa8630034adfefaf3561f3cad6ec88d323e Mon Sep 17 00:00:00 2001 From: James Hensman Date: Mon, 24 Feb 2014 08:55:18 +0000 Subject: [PATCH] some work pon EP (uninished) --- GPy/inference/latent_function_inference/ep.py | 421 +++--------------- 1 file changed, 61 insertions(+), 360 deletions(-) diff --git a/GPy/inference/latent_function_inference/ep.py b/GPy/inference/latent_function_inference/ep.py index aa106067..87c08221 100644 --- a/GPy/inference/latent_function_inference/ep.py +++ b/GPy/inference/latent_function_inference/ep.py @@ -3,390 +3,91 @@ from scipy import stats from ..util.linalg import pdinv,mdot,jitchol,chol_inv,DSYR,tdot,dtrtrs from likelihood import likelihood -class EP(likelihood): - def __init__(self,data,noise_model): - """ - Expectation Propagation - - :param data: data to model - :type data: numpy array - :param noise_model: noise distribution - :type noise_model: A GPy noise model - - """ - self.noise_model = noise_model - self.data = data - self.num_data, self.output_dim = self.data.shape - self.is_heteroscedastic = True - self.num_params = 0 - - #Initial values - Likelihood approximation parameters: - #p(y|f) = t(f|tau_tilde,v_tilde) - self.tau_tilde = np.zeros(self.num_data) - self.v_tilde = np.zeros(self.num_data) - - #initial values for the GP variables - self.Y = np.zeros((self.num_data,1)) - self.covariance_matrix = np.eye(self.num_data) - self.precision = np.ones(self.num_data)[:,None] - self.Z = 0 - self.YYT = None - self.V = self.precision * self.Y - self.VVT_factor = self.V - self.trYYT = 0. - - super(EP, self).__init__() - - def restart(self): - self.tau_tilde = np.zeros(self.num_data) - self.v_tilde = np.zeros(self.num_data) - self.Y = np.zeros((self.num_data,1)) - self.covariance_matrix = np.eye(self.num_data) - self.precision = np.ones(self.num_data)[:,None] - self.Z = 0 - self.YYT = None - self.V = self.precision * self.Y - self.VVT_factor = self.V - self.trYYT = 0. - - def predictive_values(self,mu,var,full_cov,**noise_args): - if full_cov: - raise NotImplementedError, "Cannot make correlated predictions with an EP likelihood" - return self.noise_model.predictive_values(mu,var,**noise_args) - - def log_predictive_density(self, y_test, mu_star, var_star): - """ - Calculation of the log predictive density - - .. math: - p(y_{*}|D) = p(y_{*}|f_{*})p(f_{*}|\mu_{*}\\sigma^{2}_{*}) - - :param y_test: test observations (y_{*}) - :type y_test: (Nx1) array - :param mu_star: predictive mean of gaussian p(f_{*}|mu_{*}, var_{*}) - :type mu_star: (Nx1) array - :param var_star: predictive variance of gaussian p(f_{*}|mu_{*}, var_{*}) - :type var_star: (Nx1) array - """ - return self.noise_model.log_predictive_density(y_test, mu_star, var_star) - - def _get_params(self): - #return np.zeros(0) - return self.noise_model._get_params() - - def _get_param_names(self): - #return [] - return self.noise_model._get_param_names() - - def _set_params(self,p): - #pass # TODO: the EP likelihood might want to take some parameters... - self.noise_model._set_params(p) - - def _gradients(self,partial): - #return np.zeros(0) # TODO: the EP likelihood might want to take some parameters... - return self.noise_model._gradients(partial) - - def _compute_GP_variables(self): - #Variables to be called from GP - mu_tilde = self.v_tilde/self.tau_tilde #When calling EP, this variable is used instead of Y in the GP model - sigma_sum = 1./self.tau_ + 1./self.tau_tilde - mu_diff_2 = (self.v_/self.tau_ - mu_tilde)**2 - self.Z = np.sum(np.log(self.Z_hat)) + 0.5*np.sum(np.log(sigma_sum)) + 0.5*np.sum(mu_diff_2/sigma_sum) #Normalization constant, aka Z_ep - self.Z += 0.5*self.num_data*np.log(2*np.pi) - - self.Y = mu_tilde[:,None] - self.YYT = np.dot(self.Y,self.Y.T) - self.covariance_matrix = np.diag(1./self.tau_tilde) - self.precision = self.tau_tilde[:,None] - self.V = self.precision * self.Y - self.VVT_factor = self.V - self.trYYT = np.trace(self.YYT) - - def fit_full(self, K, epsilon=1e-3,power_ep=[1.,1.]): +class EP(object): + def __init__(self, epsilon=1e-6, eta=1., delta=1.): """ The expectation-propagation algorithm. For nomenclature see Rasmussen & Williams 2006. :param epsilon: Convergence criterion, maximum squared difference allowed between mean updates to stop iterations (float) :type epsilon: float - :param power_ep: Power EP parameters - :type power_ep: list of floats - + :param eta: Power EP thing TODO: Ricardo: what, exactly? + :type eta: float64 + :param delta: Power EP thing TODO: Ricardo: what, exactly? + :type delta: float64 """ - self.epsilon = epsilon - self.eta, self.delta = power_ep + self.epsilon, self.eta, self.delta = epsilon, eta, delta + self.reset() + + def reset(self): + self.old_mutilde, self.old_vtilde = None, None + + def inference(self, kern, X, likelihood, Y, Y_metadata=None): + + K = kern.K(X) + + mu_tilde, tau_tilde = self.expectation_propagation() + + + def expectation_propagation(self, K, Y, Y_metadata, likelihood) + + num_data, data_dim = Y.shape + assert data_dim == 1, "This EP methods only works for 1D outputs" + #Initial values - Posterior distribution parameters: q(f|X,Y) = N(f|mu,Sigma) mu = np.zeros(self.num_data) Sigma = K.copy() - """ - Initial values - Cavity distribution parameters: - q_(f|mu_,sigma2_) = Product{q_i(f|mu_i,sigma2_i)} - sigma_ = 1./tau_ - mu_ = v_/tau_ - """ - self.tau_ = np.empty(self.num_data,dtype=float) - self.v_ = np.empty(self.num_data,dtype=float) - #Initial values - Marginal moments - z = np.empty(self.num_data,dtype=float) - self.Z_hat = np.empty(self.num_data,dtype=float) - phi = np.empty(self.num_data,dtype=float) - mu_hat = np.empty(self.num_data,dtype=float) - sigma2_hat = np.empty(self.num_data,dtype=float) + Z_hat = np.empty(num_data,dtype=np.float64) + mu_hat = np.empty(num_data,dtype=np.float64) + sigma2_hat = np.empty(num_data,dtype=np.float64) + + #initial values - Gaussian factors + if self.old_mutilde is None: + tau_tilde, mu_tilde, v_tilde = np.zeros((3, num_data, num_data)) + else: + assert old_mutilde.size == num_data, "data size mis-match: did you change the data? try resetting!" + mu_tilde, v_tilde = self.old_mutilde, self.old_vtilde + tau_tilde = v_tilde/mu_tilde #Approximation epsilon_np1 = self.epsilon + 1. epsilon_np2 = self.epsilon + 1. - self.iterations = 0 - self.np1 = [self.tau_tilde.copy()] - self.np2 = [self.v_tilde.copy()] - while epsilon_np1 > self.epsilon or epsilon_np2 > self.epsilon: - update_order = np.random.permutation(self.num_data) + iterations = 0 + while (epsilon_np1 > self.epsilon) or (epsilon_np2 > self.epsilon): + update_order = np.random.permutation(num_data) for i in update_order: #Cavity distribution parameters - self.tau_[i] = 1./Sigma[i,i] - self.eta*self.tau_tilde[i] - self.v_[i] = mu[i]/Sigma[i,i] - self.eta*self.v_tilde[i] + tau_cav = 1./Sigma[i,i] - self.eta*tau_tilde[i] + v_cav = mu[i]/Sigma[i,i] - self.eta*v_tilde[i] #Marginal moments - self.Z_hat[i], mu_hat[i], sigma2_hat[i] = self.noise_model.moments_match(self.data[i],self.tau_[i],self.v_[i]) + Z_hat[i], mu_hat[i], sigma2_hat[i] = likelihood.moments_match(Y[i], tau_cav, v_cav, Y_metadata=(None if Y_metadata is None else Y_metadata[i])) #Site parameters update - Delta_tau = self.delta/self.eta*(1./sigma2_hat[i] - 1./Sigma[i,i]) - Delta_v = self.delta/self.eta*(mu_hat[i]/sigma2_hat[i] - mu[i]/Sigma[i,i]) - self.tau_tilde[i] += Delta_tau - self.v_tilde[i] += Delta_v + delta_tau = self.delta/self.eta*(1./sigma2_hat[i] - 1./Sigma[i,i]) + delta_v = self.delta/self.eta*(mu_hat[i]/sigma2_hat[i] - mu[i]/Sigma[i,i]) + tau_tilde[i] += delta_tau + v_tilde[i] += delta_v #Posterior distribution parameters update - DSYR(Sigma,Sigma[:,i].copy(), -float(Delta_tau/(1.+ Delta_tau*Sigma[i,i]))) - mu = np.dot(Sigma,self.v_tilde) - self.iterations += 1 - #Sigma recomptutation with Cholesky decompositon - Sroot_tilde_K = np.sqrt(self.tau_tilde)[:,None]*K - B = np.eye(self.num_data) + np.sqrt(self.tau_tilde)[None,:]*Sroot_tilde_K + DSYR(Sigma, Sigma[:,i].copy(), -Delta_tau/(1.+ Delta_tau*Sigma[i,i])) + mu = np.dot(Sigma, v_tilde) + iterations += 1 + + #(re) compute Sigma and mu using full Cholesky decompy + tau_tilde_root = np.sqrt(tau_tilde) + Sroot_tilde_K = tau_tilde_root[:,None] * K + B = np.eye(num_data) + Sroot_tilde_K * tau_tilde_root[None,:] L = jitchol(B) - V,info = dtrtrs(L,Sroot_tilde_K,lower=1) + V, _ = dtrtrs(L, Sroot_tilde_K, lower=1) Sigma = K - np.dot(V.T,V) - mu = np.dot(Sigma,self.v_tilde) - epsilon_np1 = sum((self.tau_tilde-self.np1[-1])**2)/self.num_data - epsilon_np2 = sum((self.v_tilde-self.np2[-1])**2)/self.num_data - self.np1.append(self.tau_tilde.copy()) - self.np2.append(self.v_tilde.copy()) + mu = np.dot(Sigma,v_tilde) - return self._compute_GP_variables() + #monitor convergence + epsilon_np1 = np.mean(np.square(tau_tilde-tau_tilde_old)) + epsilon_np2 = np.mean(np.square(v_tilde-v_tilde_old)) + tau_tilde_old = tau_tilde.copy() + v_tilde_old = v_tilde.copy() - def fit_DTC(self, Kmm, Kmn, epsilon=1e-3,power_ep=[1.,1.]): - """ - The expectation-propagation algorithm with sparse pseudo-input. - For nomenclature see ... 2013. + return mu, Sigma, mu_tilde, tau_tilde - :param epsilon: Convergence criterion, maximum squared difference allowed between mean updates to stop iterations (float) - :type epsilon: float - :param power_ep: Power EP parameters - :type power_ep: list of floats - - """ - self.epsilon = epsilon - self.eta, self.delta = power_ep - - num_inducing = Kmm.shape[0] - - #TODO: this doesn't work with uncertain inputs! - - """ - Prior approximation parameters: - q(f|X) = int_{df}{N(f|KfuKuu_invu,diag(Kff-Qff)*N(u|0,Kuu)} = N(f|0,Sigma0) - Sigma0 = Qnn = Knm*Kmmi*Kmn - """ - KmnKnm = np.dot(Kmn,Kmn.T) - Lm = jitchol(Kmm) - Lmi = chol_inv(Lm) - Kmmi = np.dot(Lmi.T,Lmi) - KmmiKmn = np.dot(Kmmi,Kmn) - Qnn_diag = np.sum(Kmn*KmmiKmn,-2) - LLT0 = Kmm.copy() - - #Kmmi, Lm, Lmi, Kmm_logdet = pdinv(Kmm) - #KmnKnm = np.dot(Kmn, Kmn.T) - #KmmiKmn = np.dot(Kmmi,Kmn) - #Qnn_diag = np.sum(Kmn*KmmiKmn,-2) - #LLT0 = Kmm.copy() - - """ - Posterior approximation: q(f|y) = N(f| mu, Sigma) - Sigma = Diag + P*R.T*R*P.T + K - mu = w + P*Gamma - """ - mu = np.zeros(self.num_data) - LLT = Kmm.copy() - Sigma_diag = Qnn_diag.copy() - - """ - Initial values - Cavity distribution parameters: - q_(g|mu_,sigma2_) = Product{q_i(g|mu_i,sigma2_i)} - sigma_ = 1./tau_ - mu_ = v_/tau_ - """ - self.tau_ = np.empty(self.num_data,dtype=float) - self.v_ = np.empty(self.num_data,dtype=float) - - #Initial values - Marginal moments - z = np.empty(self.num_data,dtype=float) - self.Z_hat = np.empty(self.num_data,dtype=float) - phi = np.empty(self.num_data,dtype=float) - mu_hat = np.empty(self.num_data,dtype=float) - sigma2_hat = np.empty(self.num_data,dtype=float) - - #Approximation - epsilon_np1 = 1 - epsilon_np2 = 1 - self.iterations = 0 - np1 = [self.tau_tilde.copy()] - np2 = [self.v_tilde.copy()] - while epsilon_np1 > self.epsilon or epsilon_np2 > self.epsilon: - update_order = np.random.permutation(self.num_data) - for i in update_order: - #Cavity distribution parameters - self.tau_[i] = 1./Sigma_diag[i] - self.eta*self.tau_tilde[i] - self.v_[i] = mu[i]/Sigma_diag[i] - self.eta*self.v_tilde[i] - #Marginal moments - self.Z_hat[i], mu_hat[i], sigma2_hat[i] = self.noise_model.moments_match(self.data[i],self.tau_[i],self.v_[i]) - #Site parameters update - Delta_tau = self.delta/self.eta*(1./sigma2_hat[i] - 1./Sigma_diag[i]) - Delta_v = self.delta/self.eta*(mu_hat[i]/sigma2_hat[i] - mu[i]/Sigma_diag[i]) - self.tau_tilde[i] += Delta_tau - self.v_tilde[i] += Delta_v - #Posterior distribution parameters update - DSYR(LLT,Kmn[:,i].copy(),Delta_tau) #LLT = LLT + np.outer(Kmn[:,i],Kmn[:,i])*Delta_tau - L = jitchol(LLT) - #cholUpdate(L,Kmn[:,i]*np.sqrt(Delta_tau)) - V,info = dtrtrs(L,Kmn,lower=1) - Sigma_diag = np.sum(V*V,-2) - si = np.sum(V.T*V[:,i],-1) - mu += (Delta_v-Delta_tau*mu[i])*si - self.iterations += 1 - #Sigma recomputation with Cholesky decompositon - LLT = LLT0 + np.dot(Kmn*self.tau_tilde[None,:],Kmn.T) - L = jitchol(LLT) - V,info = dtrtrs(L,Kmn,lower=1) - V2,info = dtrtrs(L.T,V,lower=0) - Sigma_diag = np.sum(V*V,-2) - Knmv_tilde = np.dot(Kmn,self.v_tilde) - mu = np.dot(V2.T,Knmv_tilde) - epsilon_np1 = sum((self.tau_tilde-np1[-1])**2)/self.num_data - epsilon_np2 = sum((self.v_tilde-np2[-1])**2)/self.num_data - np1.append(self.tau_tilde.copy()) - np2.append(self.v_tilde.copy()) - - self._compute_GP_variables() - - def fit_FITC(self, Kmm, Kmn, Knn_diag, epsilon=1e-3,power_ep=[1.,1.]): - """ - The expectation-propagation algorithm with sparse pseudo-input. - For nomenclature see Naish-Guzman and Holden, 2008. - - :param epsilon: Convergence criterion, maximum squared difference allowed between mean updates to stop iterations (float) - :type epsilon: float - :param power_ep: Power EP parameters - :type power_ep: list of floats - """ - self.epsilon = epsilon - self.eta, self.delta = power_ep - - num_inducing = Kmm.shape[0] - - """ - Prior approximation parameters: - q(f|X) = int_{df}{N(f|KfuKuu_invu,diag(Kff-Qff)*N(u|0,Kuu)} = N(f|0,Sigma0) - Sigma0 = diag(Knn-Qnn) + Qnn, Qnn = Knm*Kmmi*Kmn - """ - Lm = jitchol(Kmm) - Lmi = chol_inv(Lm) - Kmmi = np.dot(Lmi.T,Lmi) - P0 = Kmn.T - KmnKnm = np.dot(P0.T, P0) - KmmiKmn = np.dot(Kmmi,P0.T) - Qnn_diag = np.sum(P0.T*KmmiKmn,-2) - Diag0 = Knn_diag - Qnn_diag - R0 = jitchol(Kmmi).T - - """ - Posterior approximation: q(f|y) = N(f| mu, Sigma) - Sigma = Diag + P*R.T*R*P.T + K - mu = w + P*Gamma - """ - self.w = np.zeros(self.num_data) - self.Gamma = np.zeros(num_inducing) - mu = np.zeros(self.num_data) - P = P0.copy() - R = R0.copy() - Diag = Diag0.copy() - Sigma_diag = Knn_diag - RPT0 = np.dot(R0,P0.T) - - """ - Initial values - Cavity distribution parameters: - q_(g|mu_,sigma2_) = Product{q_i(g|mu_i,sigma2_i)} - sigma_ = 1./tau_ - mu_ = v_/tau_ - """ - self.tau_ = np.empty(self.num_data,dtype=float) - self.v_ = np.empty(self.num_data,dtype=float) - - #Initial values - Marginal moments - z = np.empty(self.num_data,dtype=float) - self.Z_hat = np.empty(self.num_data,dtype=float) - phi = np.empty(self.num_data,dtype=float) - mu_hat = np.empty(self.num_data,dtype=float) - sigma2_hat = np.empty(self.num_data,dtype=float) - - #Approximation - epsilon_np1 = 1 - epsilon_np2 = 1 - self.iterations = 0 - self.np1 = [self.tau_tilde.copy()] - self.np2 = [self.v_tilde.copy()] - while epsilon_np1 > self.epsilon or epsilon_np2 > self.epsilon: - update_order = np.random.permutation(self.num_data) - for i in update_order: - #Cavity distribution parameters - self.tau_[i] = 1./Sigma_diag[i] - self.eta*self.tau_tilde[i] - self.v_[i] = mu[i]/Sigma_diag[i] - self.eta*self.v_tilde[i] - #Marginal moments - self.Z_hat[i], mu_hat[i], sigma2_hat[i] = self.noise_model.moments_match(self.data[i],self.tau_[i],self.v_[i]) - #Site parameters update - Delta_tau = self.delta/self.eta*(1./sigma2_hat[i] - 1./Sigma_diag[i]) - Delta_v = self.delta/self.eta*(mu_hat[i]/sigma2_hat[i] - mu[i]/Sigma_diag[i]) - self.tau_tilde[i] += Delta_tau - self.v_tilde[i] += Delta_v - #Posterior distribution parameters update - dtd1 = Delta_tau*Diag[i] + 1. - dii = Diag[i] - Diag[i] = dii - (Delta_tau * dii**2.)/dtd1 - pi_ = P[i,:].reshape(1,num_inducing) - P[i,:] = pi_ - (Delta_tau*dii)/dtd1 * pi_ - Rp_i = np.dot(R,pi_.T) - RTR = np.dot(R.T,np.dot(np.eye(num_inducing) - Delta_tau/(1.+Delta_tau*Sigma_diag[i]) * np.dot(Rp_i,Rp_i.T),R)) - R = jitchol(RTR).T - self.w[i] += (Delta_v - Delta_tau*self.w[i])*dii/dtd1 - self.Gamma += (Delta_v - Delta_tau*mu[i])*np.dot(RTR,P[i,:].T) - RPT = np.dot(R,P.T) - Sigma_diag = Diag + np.sum(RPT.T*RPT.T,-1) - mu = self.w + np.dot(P,self.Gamma) - self.iterations += 1 - #Sigma recomptutation with Cholesky decompositon - Iplus_Dprod_i = 1./(1.+ Diag0 * self.tau_tilde) - Diag = Diag0 * Iplus_Dprod_i - P = Iplus_Dprod_i[:,None] * P0 - safe_diag = np.where(Diag0 < self.tau_tilde, self.tau_tilde/(1.+Diag0*self.tau_tilde), (1. - Iplus_Dprod_i)/Diag0) - L = jitchol(np.eye(num_inducing) + np.dot(RPT0,safe_diag[:,None]*RPT0.T)) - R,info = dtrtrs(L,R0,lower=1) - RPT = np.dot(R,P.T) - Sigma_diag = Diag + np.sum(RPT.T*RPT.T,-1) - self.w = Diag * self.v_tilde - self.Gamma = np.dot(R.T, np.dot(RPT,self.v_tilde)) - mu = self.w + np.dot(P,self.Gamma) - epsilon_np1 = sum((self.tau_tilde-self.np1[-1])**2)/self.num_data - epsilon_np2 = sum((self.v_tilde-self.np2[-1])**2)/self.num_data - self.np1.append(self.tau_tilde.copy()) - self.np2.append(self.v_tilde.copy()) - - return self._compute_GP_variables()