From eee1fba8adee3f2c4601f9888c24cb8f5dcb188a Mon Sep 17 00:00:00 2001 From: Thang Bui Date: Wed, 19 Oct 2016 15:56:10 +0100 Subject: [PATCH 1/4] Added pep.py -- Sparse Gaussian processes using Power Expectation Propagation This allows interpolation between FITC (EP or alpha = 1), and Titsias's variational (VarDTC, VFE when alpha = 0). Reference: A Unifying Framework for Sparse Gaussian Process Approximation using Power Expectation Propagation https://arxiv.org/abs/1605.07066 --- .../latent_function_inference/pep.py | 93 +++++++++++++++++++ 1 file changed, 93 insertions(+) create mode 100644 GPy/inference/latent_function_inference/pep.py diff --git a/GPy/inference/latent_function_inference/pep.py b/GPy/inference/latent_function_inference/pep.py new file mode 100644 index 00000000..79706292 --- /dev/null +++ b/GPy/inference/latent_function_inference/pep.py @@ -0,0 +1,93 @@ +from .posterior import Posterior +from ...util.linalg import jitchol, tdot, dtrtrs, dtrtri, pdinv +from ...util import diag +import numpy as np +from . import LatentFunctionInference +log_2_pi = np.log(2*np.pi) + +class PEP(LatentFunctionInference): + ''' + Sparse Gaussian processes using Power-Expectation Propagation + for regression: alpha \approx 0 gives VarDTC and alpha = 1 gives FITC + + Reference: A Unifying Framework for Sparse Gaussian Process Approximation using + Power Expectation Propagation, https://arxiv.org/abs/1605.07066 + + ''' + const_jitter = 1e-6 + + def __init__(self, alpha): + super(PEP, self).__init__() + self.alpha = alpha + + def inference(self, kern, X, Z, likelihood, Y, mean_function=None, Y_metadata=None): + assert mean_function is None, "inference with a mean function not implemented" + + num_inducing, _ = Z.shape + num_data, output_dim = Y.shape + + #make sure the noise is not hetero + sigma_n = likelihood.gaussian_variance(Y_metadata) + if sigma_n.size >1: + raise NotImplementedError("no hetero noise with this implementation of PEP") + + Kmm = kern.K(Z) + Knn = kern.Kdiag(X) + Knm = kern.K(X, Z) + U = Knm + + #factor Kmm + diag.add(Kmm, self.const_jitter) + Kmmi, L, Li, _ = pdinv(Kmm) + + #compute beta_star, the effective noise precision + LiUT = np.dot(Li, U.T) + sigma_star = sigma_n + self.alpha * (Knn - np.sum(np.square(LiUT),0)) + beta_star = 1./sigma_star + + # Compute and factor A + A = tdot(LiUT*np.sqrt(beta_star)) + np.eye(num_inducing) + LA = jitchol(A) + + # back substitute to get b, P, v + URiy = np.dot(U.T*beta_star,Y) + tmp, _ = dtrtrs(L, URiy, lower=1) + b, _ = dtrtrs(LA, tmp, lower=1) + tmp, _ = dtrtrs(LA, b, lower=1, trans=1) + v, _ = dtrtrs(L, tmp, lower=1, trans=1) + tmp, _ = dtrtrs(LA, Li, lower=1, trans=0) + P = tdot(tmp.T) + + alpha_const_term = (1.0-self.alpha) / self.alpha + + #compute log marginal + log_marginal = -0.5*num_data*output_dim*np.log(2*np.pi) + \ + -np.sum(np.log(np.diag(LA)))*output_dim + \ + 0.5*output_dim*(1+alpha_const_term)*np.sum(np.log(beta_star)) + \ + -0.5*np.sum(np.square(Y.T*np.sqrt(beta_star))) + \ + 0.5*np.sum(np.square(b)) + 0.5*alpha_const_term*num_data*np.log(sigma_n) + #compute dL_dR + Uv = np.dot(U, v) + dL_dR = 0.5*(np.sum(U*np.dot(U,P), 1) - (1.0+alpha_const_term)/beta_star + np.sum(np.square(Y), 1) - 2.*np.sum(Uv*Y, 1) \ + + np.sum(np.square(Uv), 1))*beta_star**2 + + # Compute dL_dKmm + vvT_P = tdot(v.reshape(-1,1)) + P + dL_dK = 0.5*(Kmmi - vvT_P) + KiU = np.dot(Kmmi, U.T) + dL_dK += self.alpha * np.dot(KiU*dL_dR, KiU.T) + + # Compute dL_dU + vY = np.dot(v.reshape(-1,1),Y.T) + dL_dU = vY - np.dot(vvT_P, U.T) + dL_dU *= beta_star + dL_dU -= self.alpha * 2.*KiU*dL_dR + + dL_dthetaL = likelihood.exact_inference_gradients(dL_dR) + dL_dthetaL += 0.5*alpha_const_term*num_data / sigma_n + grad_dict = {'dL_dKmm': dL_dK, 'dL_dKdiag':dL_dR * self.alpha, 'dL_dKnm':dL_dU.T, 'dL_dthetaL':dL_dthetaL} + + #construct a posterior object + post = Posterior(woodbury_inv=Kmmi-P, woodbury_vector=v, K=Kmm, mean=None, cov=None, K_chol=L) + + return post, log_marginal, grad_dict From df7c7539f924af0a21b909b55a1bf1d3899d7790 Mon Sep 17 00:00:00 2001 From: Thang Bui Date: Thu, 27 Oct 2016 14:39:32 +0100 Subject: [PATCH 2/4] added tests --- .../latent_function_inference/__init__.py | 1 + GPy/testing/pep.py | 92 +++++++++++++++++++ 2 files changed, 93 insertions(+) create mode 100644 GPy/testing/pep.py diff --git a/GPy/inference/latent_function_inference/__init__.py b/GPy/inference/latent_function_inference/__init__.py index 90fbf0f1..3938a6a4 100644 --- a/GPy/inference/latent_function_inference/__init__.py +++ b/GPy/inference/latent_function_inference/__init__.py @@ -67,6 +67,7 @@ from GPy.inference.latent_function_inference.var_dtc import VarDTC from .expectation_propagation import EP, EPDTC from .dtc import DTC from .fitc import FITC +from .pep import PEP from .var_dtc_parallel import VarDTC_minibatch from .var_gauss import VarGauss from .gaussian_grid_inference import GaussianGridInference diff --git a/GPy/testing/pep.py b/GPy/testing/pep.py new file mode 100644 index 00000000..0d159d5d --- /dev/null +++ b/GPy/testing/pep.py @@ -0,0 +1,92 @@ +# Copyright (c) 2014, James Hensman, 2016, Thang Bui +# Licensed under the BSD 3-clause license (see LICENSE.txt) + +import unittest +import numpy as np +import GPy + +class PEPgradienttest(unittest.TestCase): + def setUp(self): + ###################################### + # # 1 dimensional example + + N = 20 + # sample inputs and outputs + self.X1D = np.random.uniform(-3., 3., (N, 1)) + self.Y1D = np.sin(self.X1D) + np.random.randn(N, 1) * 0.05 + + ###################################### + # # 2 dimensional example + + # sample inputs and outputs + self.X2D = np.random.uniform(-3., 3., (N, 2)) + self.Y2D = np.sin(self.X2D[:, 0:1]) * np.sin(self.X2D[:, 1:2]) + np.random.randn(N, 1) * 0.05 + + ####################################### + # # more datapoints, check in alpha limits, the log marginal likelihood + # # is consistent with FITC and VFE/Var_DTC + M = 5 + self.X1 = np.c_[np.linspace(-1., 1., N)] + self.Y1 = np.sin(self.X1) + np.random.randn(N, 1) * 0.05 + self.kernel = GPy.kern.RBF(input_dim=1, lengthscale=0.5, variance=1) + self.Z = np.random.uniform(-1, 1, (M, 1)) + self.lik_noise_var = 0.01 + + def test_pep_1d_gradients(self): + m = GPy.models.SparseGPRegression(self.X1D, self.Y1D) + m.inference_method = GPy.inference.latent_function_inference.PEP(alpha=np.random.rand()) + self.assertTrue(m.checkgrad()) + + def test_pep_2d_gradients(self): + m = GPy.models.SparseGPRegression(self.X2D, self.Y2D) + m.inference_method = GPy.inference.latent_function_inference.PEP(alpha=np.random.rand()) + self.assertTrue(m.checkgrad()) + + def test_pep_vfe_consistency(self): + vfe_model = GPy.models.SparseGPRegression( + self.X1, + self.Y1, + kernel=self.kernel, + Z=self.Z + ) + vfe_model.inference_method = GPy.inference.latent_function_inference.VarDTC() + vfe_model.Gaussian_noise.variance = self.lik_noise_var + vfe_lml = vfe_model.log_likelihood() + + pep_model = GPy.models.SparseGPRegression( + self.X1, + self.Y1, + kernel=self.kernel, + Z=self.Z + ) + pep_model.inference_method = GPy.inference.latent_function_inference.PEP(alpha=1e-5) + pep_model.Gaussian_noise.variance = self.lik_noise_var + pep_lml = pep_model.log_likelihood() + + self.assertAlmostEqual(vfe_lml[0, 0], pep_lml[0], delta=abs(0.01*pep_lml[0])) + + def test_pep_fitc_consistency(self): + fitc_model = GPy.models.SparseGPRegression( + self.X1D, + self.Y1D, + kernel=self.kernel, + Z=self.Z + ) + fitc_model.inference_method = GPy.inference.latent_function_inference.FITC() + fitc_model.Gaussian_noise.variance = self.lik_noise_var + fitc_lml = fitc_model.log_likelihood() + + pep_model = GPy.models.SparseGPRegression( + self.X1D, + self.Y1D, + kernel=self.kernel, + Z=self.Z + ) + pep_model.inference_method = GPy.inference.latent_function_inference.PEP(alpha=1) + pep_model.Gaussian_noise.variance = self.lik_noise_var + pep_lml = pep_model.log_likelihood() + + self.assertAlmostEqual(fitc_lml, pep_lml[0], delta=abs(0.001*pep_lml[0])) + + + From 323204541e04c5c3ca981188a2d04f4a6a1721ef Mon Sep 17 00:00:00 2001 From: Thang Bui Date: Tue, 1 Nov 2016 16:48:45 +0000 Subject: [PATCH 3/4] fixed seed in pep test script #448 --- GPy/testing/pep.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/GPy/testing/pep.py b/GPy/testing/pep.py index 0d159d5d..2aa6a784 100644 --- a/GPy/testing/pep.py +++ b/GPy/testing/pep.py @@ -9,6 +9,7 @@ class PEPgradienttest(unittest.TestCase): def setUp(self): ###################################### # # 1 dimensional example + np.random.seed(10) N = 20 # sample inputs and outputs @@ -26,6 +27,7 @@ class PEPgradienttest(unittest.TestCase): # # more datapoints, check in alpha limits, the log marginal likelihood # # is consistent with FITC and VFE/Var_DTC M = 5 + np.random.seed(42) self.X1 = np.c_[np.linspace(-1., 1., N)] self.Y1 = np.sin(self.X1) + np.random.randn(N, 1) * 0.05 self.kernel = GPy.kern.RBF(input_dim=1, lengthscale=0.5, variance=1) From 7b384ae25878171815ccd164dcbddea8b0ec0f2c Mon Sep 17 00:00:00 2001 From: Thang Bui Date: Mon, 7 Nov 2016 09:54:50 +0000 Subject: [PATCH 4/4] renamed pep test scripts --- GPy/testing/{pep.py => pep_tests.py} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename GPy/testing/{pep.py => pep_tests.py} (100%) diff --git a/GPy/testing/pep.py b/GPy/testing/pep_tests.py similarity index 100% rename from GPy/testing/pep.py rename to GPy/testing/pep_tests.py