From 167e9c538dc6862a2581dbbe9c8e02d1765052b1 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Fri, 30 Jan 2015 09:43:02 +0000 Subject: [PATCH 1/5] [var dtc] added code for additional covariates, not affecting normal procedures --- GPy/models/bayesian_gplvm_minibatch.py | 4 ++-- GPy/models/sparse_gp_minibatch.py | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/GPy/models/bayesian_gplvm_minibatch.py b/GPy/models/bayesian_gplvm_minibatch.py index 64aed246..71f69eb2 100644 --- a/GPy/models/bayesian_gplvm_minibatch.py +++ b/GPy/models/bayesian_gplvm_minibatch.py @@ -83,8 +83,8 @@ class BayesianGPLVMMiniBatch(SparseGPMiniBatch): """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) + def _inner_parameters_changed(self, kern, X, Z, likelihood, Y, Y_metadata, Lm=None, dL_dKmm=None, subset_indices=None, **kw): + 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, **kw) if self.has_uncertain_inputs(): current_values['meangrad'], current_values['vargrad'] = self.kern.gradients_qX_expectations( diff --git a/GPy/models/sparse_gp_minibatch.py b/GPy/models/sparse_gp_minibatch.py index ec2e28f5..8925d4d7 100644 --- a/GPy/models/sparse_gp_minibatch.py +++ b/GPy/models/sparse_gp_minibatch.py @@ -97,7 +97,7 @@ Created on 3 Nov 2014 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): + def _inner_parameters_changed(self, kern, X, Z, likelihood, Y, Y_metadata, Lm=None, dL_dKmm=None, subset_indices=None, **kwargs): """ This is the standard part, which usually belongs in parameters_changed. @@ -117,7 +117,7 @@ Created on 3 Nov 2014 algorithm. """ try: - posterior, log_marginal_likelihood, grad_dict = self.inference_method.inference(kern, X, Z, likelihood, Y, Y_metadata, Lm=Lm, dL_dKmm=None) + posterior, log_marginal_likelihood, grad_dict = self.inference_method.inference(kern, X, Z, likelihood, Y, Y_metadata, Lm=Lm, dL_dKmm=None, **kwargs) except: posterior, log_marginal_likelihood, grad_dict = self.inference_method.inference(kern, X, Z, likelihood, Y, Y_metadata) current_values = {} From 29d153e185bb1a42f27dff2030020d75dc9026bb Mon Sep 17 00:00:00 2001 From: Alan Saul Date: Mon, 9 Feb 2015 19:35:46 +0000 Subject: [PATCH 2/5] Bug in linalg jitchol!!! --- GPy/testing/linalg_test.py | 35 +++++++++++++++++++++++++++++++++++ GPy/util/linalg.py | 12 ++++++------ 2 files changed, 41 insertions(+), 6 deletions(-) create mode 100644 GPy/testing/linalg_test.py diff --git a/GPy/testing/linalg_test.py b/GPy/testing/linalg_test.py new file mode 100644 index 00000000..b734f6af --- /dev/null +++ b/GPy/testing/linalg_test.py @@ -0,0 +1,35 @@ +import numpy as np +import scipy as sp +from ..util.linalg import jitchol + +class LinalgTests(np.testing.TestCase): + def setUp(self): + #Create PD matrix + A = np.random.randn(20,100) + self.A = A.dot(A.T) + #compute Eigdecomp + vals, vectors = np.linalg.eig(self.A) + #Set smallest eigenval to be negative with 5 rounds worth of jitter + vals[vals.argmin()] = 0 + default_jitter = 1e-6*np.mean(vals) + vals[vals.argmin()] = -default_jitter*(10**3.5) + self.A_corrupt = (vectors * vals).dot(vectors.T) + + def test_jitchol_success(self): + """ + Expect 5 rounds of jitter to be added and for the recovered matrix to be + identical to the corrupted matrix apart from the jitter added to the diagonal + """ + L = jitchol(self.A_corrupt, maxtries=5) + A_new = L.dot(L.T) + diff = A_new - self.A_corrupt + np.testing.assert_allclose(diff, np.eye(A_new.shape[0])*np.diag(diff).mean(), atol=1e-13) + + def test_jitchol_failure(self): + try: + """ Expecting an exception to be thrown as we expect it to require + 5 rounds of jitter to be added to enforce PDness""" + jitchol(self.A_corrupt, maxtries=4) + return False + except sp.linalg.LinAlgError: + return True diff --git a/GPy/util/linalg.py b/GPy/util/linalg.py index dffd438a..2c02357c 100644 --- a/GPy/util/linalg.py +++ b/GPy/util/linalg.py @@ -82,6 +82,7 @@ def force_F_ordered(A): # return jitchol(A+np.eye(A.shape[0])*jitter, maxtries-1) + def jitchol(A, maxtries=5): A = np.ascontiguousarray(A) L, info = lapack.dpotrf(A, lower=1) @@ -92,13 +93,16 @@ def jitchol(A, maxtries=5): if np.any(diagA <= 0.): raise linalg.LinAlgError, "not pd: non-positive diagonal elements" jitter = diagA.mean() * 1e-6 - while maxtries > 0 and np.isfinite(jitter): + num_tries = 0 + while num_tries < maxtries and np.isfinite(jitter): try: + print jitter L = linalg.cholesky(A + np.eye(A.shape[0]) * jitter, lower=True) + return L except: jitter *= 10 finally: - maxtries -= 1 + num_tries += 1 raise linalg.LinAlgError, "not positive definite, even with jitter." import traceback try: raise @@ -108,10 +112,6 @@ def jitchol(A, maxtries=5): import ipdb;ipdb.set_trace() return L - - - - # def dtrtri(L, lower=1): # """ # Wrapper for lapack dtrtri function From f69019238406c427843abdc5977d2dd3b21b4a20 Mon Sep 17 00:00:00 2001 From: Alan Saul Date: Tue, 10 Feb 2015 11:52:40 +0000 Subject: [PATCH 3/5] Added logging for jitter so we know how much has been added and how many tries have been taken --- GPy/testing/linalg_test.py | 6 ++++-- GPy/util/linalg.py | 17 ++++++----------- 2 files changed, 10 insertions(+), 13 deletions(-) diff --git a/GPy/testing/linalg_test.py b/GPy/testing/linalg_test.py index b734f6af..8e103795 100644 --- a/GPy/testing/linalg_test.py +++ b/GPy/testing/linalg_test.py @@ -27,8 +27,10 @@ class LinalgTests(np.testing.TestCase): def test_jitchol_failure(self): try: - """ Expecting an exception to be thrown as we expect it to require - 5 rounds of jitter to be added to enforce PDness""" + """ + Expecting an exception to be thrown as we expect it to require + 5 rounds of jitter to be added to enforce PDness + """ jitchol(self.A_corrupt, maxtries=4) return False except sp.linalg.LinAlgError: diff --git a/GPy/util/linalg.py b/GPy/util/linalg.py index 2c02357c..b148f2f4 100644 --- a/GPy/util/linalg.py +++ b/GPy/util/linalg.py @@ -93,24 +93,19 @@ def jitchol(A, maxtries=5): if np.any(diagA <= 0.): raise linalg.LinAlgError, "not pd: non-positive diagonal elements" jitter = diagA.mean() * 1e-6 - num_tries = 0 - while num_tries < maxtries and np.isfinite(jitter): + num_tries = 1 + while num_tries <= maxtries and np.isfinite(jitter): try: - print jitter L = linalg.cholesky(A + np.eye(A.shape[0]) * jitter, lower=True) + logging.warning('Added {} rounds of jitter, jitter of {:.10e}\n'.format(num_tries, jitter)) return L except: jitter *= 10 - finally: num_tries += 1 - raise linalg.LinAlgError, "not positive definite, even with jitter." import traceback - try: raise - except: - logging.warning('\n'.join(['Added jitter of {:.10e}'.format(jitter), - ' in '+traceback.format_list(traceback.extract_stack(limit=2)[-2:-1])[0][2:]])) - import ipdb;ipdb.set_trace() - return L + logging.warning('\n'.join(['Added {} rounds of jitter, jitter of {:.10e}'.format(num_tries-1, jitter), + ' in '+traceback.format_list(traceback.extract_stack(limit=2)[-2:-1])[0][2:]])) + raise linalg.LinAlgError, "not positive definite, even with jitter." # def dtrtri(L, lower=1): # """ From b499a870fbce243d7989ab8710c8fa1d83af00ea Mon Sep 17 00:00:00 2001 From: James Hensman Date: Tue, 10 Feb 2015 14:23:03 +0000 Subject: [PATCH 4/5] fixed a plotting bug for sliced plots --- GPy/likelihoods/likelihood.py | 2 +- GPy/plotting/matplot_dep/models_plots.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/GPy/likelihoods/likelihood.py b/GPy/likelihoods/likelihood.py index 87b7315e..790c6ba4 100644 --- a/GPy/likelihoods/likelihood.py +++ b/GPy/likelihoods/likelihood.py @@ -178,7 +178,7 @@ class Likelihood(Parameterized): stop dF_dtheta = None # Not yet implemented - return F.reshape(*shape), dF_dm.reshape(*shape), dF_dv.reshape(*shape), None + return F.reshape(*shape), dF_dm.reshape(*shape), dF_dv.reshape(*shape), dF_dtheta def predictive_mean(self, mu, variance, Y_metadata=None): """ diff --git a/GPy/plotting/matplot_dep/models_plots.py b/GPy/plotting/matplot_dep/models_plots.py index ed024c0a..d2d5a8e2 100644 --- a/GPy/plotting/matplot_dep/models_plots.py +++ b/GPy/plotting/matplot_dep/models_plots.py @@ -172,7 +172,7 @@ def plot_fit(model, plot_limits=None, which_data_rows='all', if hasattr(model,"Z"): #Zu = model.Z[:,free_dims] * model._Xscale[:,free_dims] + model._Xoffset[:,free_dims] Zu = Z[:,free_dims] - plots['inducing_inputs'] = ax.plot(Zu[:,free_dims[0]], Zu[:,free_dims[1]], 'wo') + plots['inducing_inputs'] = ax.plot(Zu[:,0], Zu[:,1], 'wo') else: raise NotImplementedError, "Cannot define a frame with more than two input dimensions" From 06d7e690f371277ba40c89a7785f9dcb4db69c81 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Tue, 10 Feb 2015 15:54:59 +0000 Subject: [PATCH 5/5] minor weave/numpy bug in coregionalize --- GPy/kern/_src/coregionalize.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/GPy/kern/_src/coregionalize.py b/GPy/kern/_src/coregionalize.py index 3fcf1c98..cbfe1ccb 100644 --- a/GPy/kern/_src/coregionalize.py +++ b/GPy/kern/_src/coregionalize.py @@ -127,7 +127,7 @@ class Coregionalize(Kern): config.set('weave', 'working', 'False') dL_dK_small = self._gradient_reduce_weave(dL_dK, index, index2) else: - dL_dK_small = self._gradient_reduce_weave(dL_dK, index, index2) + dL_dK_small = self._gradient_reduce_numpy(dL_dK, index, index2)