From 70e4d8aee5070b3b6fdf7e4739b42178a2d34f28 Mon Sep 17 00:00:00 2001 From: Joachim van der Herten Date: Fri, 5 May 2017 10:59:15 +0200 Subject: [PATCH 1/5] Adjusted model initial hyperparameters in the tests. Also added functionality to restart from the initially supplied hyperparameters. --- GPflowOpt/acquisition.py | 14 +++++--- testing/test_acquisition.py | 68 +++++++++++++++++++++---------------- testing/test_optimizers.py | 2 +- 3 files changed, 49 insertions(+), 35 deletions(-) diff --git a/GPflowOpt/acquisition.py b/GPflowOpt/acquisition.py index c24e553..e8f6c21 100644 --- a/GPflowOpt/acquisition.py +++ b/GPflowOpt/acquisition.py @@ -40,14 +40,18 @@ class Acquisition(Parameterized): def __init__(self, models=[]): super(Acquisition, self).__init__() self.models = ParamList(np.atleast_1d(models).tolist()) + self._default_params = map(lambda m: m.get_free_state(), self.models) self._optimize_all() def _optimize_all(self): - for model in self.models: - # If likelihood variance is close to zero, updating data may result in non-invertible K - # Increase likelihood variance a bit. - model.likelihood.variance = 4.0 - model.optimize() + for model, hypers in zip(self.models, self._default_params): + try: + model.optimize() + except: + # After data update, the starting point for the hypers may result in non-invertible K + # Reset the hypers to the values when the acquisition function was initialized. + model.set_state(hypers) + model.optimize() def _build_acquisition_wrapper(self, Xcand, gradients=True): acq = self.build_acquisition(Xcand) diff --git a/testing/test_acquisition.py b/testing/test_acquisition.py index 651fee2..84af083 100644 --- a/testing/test_acquisition.py +++ b/testing/test_acquisition.py @@ -22,6 +22,20 @@ class _TestAcquisition(object): def domain(self): return np.sum([GPflowOpt.domain.ContinuousParameter("x{0}".format(i), -1, 1) for i in range(1, 3)]) + def create_parabola_model(self, design=None): + if design is None: + design = GPflowOpt.design.FactorialDesign(4, self.domain) + X, Y = design.generate(), parabola2d(design.generate()) + m = GPflow.gpr.GPR(X, Y, GPflow.kernels.RBF(2, ARD=True, lengthscales=X.std(axis=0))) + return m + + def create_plane_model(self, design=None): + if design is None: + design = GPflowOpt.design.FactorialDesign(5, self.domain) + X, Y = design.generate(), plane(design.generate()) + m = GPflow.gpr.GPR(X, Y, GPflow.kernels.RBF(2, ARD=True, lengthscales=X.std(axis=0))) + return m + def setUp(self): self.acquisition = None @@ -73,9 +87,7 @@ def test_data_update(self): class TestExpectedImprovement(_TestAcquisition, unittest.TestCase): def setUp(self): super(TestExpectedImprovement, self).setUp() - design = GPflowOpt.design.FactorialDesign(4, self.domain) - X, Y = design.generate(), parabola2d(design.generate()) - self.model = GPflow.gpr.GPR(X, Y, GPflow.kernels.RBF(2, ARD=True)) + self.model = self.create_parabola_model() self.acquisition = GPflowOpt.acquisition.ExpectedImprovement(self.model) def test_object_integrity(self): @@ -109,9 +121,7 @@ def test_EI_validity(self): class TestProbabilityOfFeasibility(_TestAcquisition, unittest.TestCase): def setUp(self): super(TestProbabilityOfFeasibility, self).setUp() - design = GPflowOpt.design.FactorialDesign(5, self.domain) - X, Y = design.generate(), plane(design.generate()) - self.model = GPflow.gpr.GPR(X, Y, GPflow.kernels.RBF(2, ARD=True)) + self.model = self.create_plane_model() self.acquisition = GPflowOpt.acquisition.ProbabilityOfFeasibility(self.model) def test_object_integrity(self): @@ -128,7 +138,6 @@ def test_PoF_validity(self): class _TestAcquisitionBinaryOperator(_TestAcquisition): - def test_object_integrity(self): self.assertTrue(isinstance(self.acquisition.lhs, GPflowOpt.acquisition.Acquisition), msg="Left hand operand should be an acquisition object") @@ -150,22 +159,23 @@ def test_data(self): class TestAcquisitionSum(_TestAcquisitionBinaryOperator, unittest.TestCase): def setUp(self): super(TestAcquisitionSum, self).setUp() - design = GPflowOpt.design.FactorialDesign(4, self.domain) - X, Y = design.generate(), parabola2d(design.generate()) - self.model = GPflow.gpr.GPR(X, Y, GPflow.kernels.RBF(2, ARD=True)) - self.acquisition = GPflowOpt.acquisition.AcquisitionSum(GPflowOpt.acquisition.ExpectedImprovement(self.model), - GPflowOpt.acquisition.ExpectedImprovement(self.model)) + self.models = [self.create_parabola_model(), self.create_parabola_model()] + self.acquisition = GPflowOpt.acquisition.AcquisitionSum( + GPflowOpt.acquisition.ExpectedImprovement(self.models[0]), + GPflowOpt.acquisition.ExpectedImprovement(self.models[1])) def test_sum_validity(self): - single_ei = GPflowOpt.acquisition.ExpectedImprovement(self.model) design = GPflowOpt.design.FactorialDesign(4, self.domain) + m = self.create_parabola_model(design) + single_ei = GPflowOpt.acquisition.ExpectedImprovement(m) p1 = self.acquisition.evaluate(design.generate()) p2 = single_ei.evaluate(design.generate()) np.testing.assert_allclose(p2, p1 / 2, rtol=1e-3, err_msg="The sum of 2 EI should be the double of only EI") def test_generating_operator(self): - joint = GPflowOpt.acquisition.ExpectedImprovement(self.model) + \ - GPflowOpt.acquisition.ExpectedImprovement(self.model) + design = GPflowOpt.design.FactorialDesign(4, self.domain) + joint = GPflowOpt.acquisition.ExpectedImprovement(self.create_parabola_model(design)) + \ + GPflowOpt.acquisition.ExpectedImprovement(self.create_parabola_model(design)) self.assertTrue(isinstance(joint, GPflowOpt.acquisition.AcquisitionSum)) def test_indices(self): @@ -178,24 +188,23 @@ def test_indices(self): class TestAcquisitionProduct(_TestAcquisitionBinaryOperator, unittest.TestCase): def setUp(self): super(TestAcquisitionProduct, self).setUp() - design = GPflowOpt.design.FactorialDesign(4, self.domain) - X, Y = design.generate(), parabola2d(design.generate()) - self.model = GPflow.gpr.GPR(X, Y, GPflow.kernels.RBF(2, ARD=True)) + self.models = [self.create_parabola_model(), self.create_parabola_model()] self.acquisition = GPflowOpt.acquisition.AcquisitionProduct( - GPflowOpt.acquisition.ExpectedImprovement(self.model), - GPflowOpt.acquisition.ExpectedImprovement(self.model)) + GPflowOpt.acquisition.ExpectedImprovement(self.models[0]), + GPflowOpt.acquisition.ExpectedImprovement(self.models[1])) def test_product_validity(self): - single_ei = GPflowOpt.acquisition.ExpectedImprovement(self.model) design = GPflowOpt.design.FactorialDesign(4, self.domain) + m = self.create_parabola_model(design) + single_ei = GPflowOpt.acquisition.ExpectedImprovement(m) p1 = self.acquisition.evaluate(design.generate()) p2 = single_ei.evaluate(design.generate()) np.testing.assert_allclose(p2, np.sqrt(p1), rtol=1e-3, err_msg="The product of 2 EI should be the square of one EI") def test_generating_operator(self): - joint = GPflowOpt.acquisition.ExpectedImprovement(self.model) * \ - GPflowOpt.acquisition.ExpectedImprovement(self.model) + joint = GPflowOpt.acquisition.ExpectedImprovement(self.create_parabola_model()) * \ + GPflowOpt.acquisition.ExpectedImprovement(self.create_parabola_model()) self.assertTrue(isinstance(joint, GPflowOpt.acquisition.AcquisitionProduct)) def test_indices(self): @@ -206,7 +215,6 @@ def test_indices(self): class TestJointAcquisition(unittest.TestCase): - @property def domain(self): return np.sum([GPflowOpt.domain.ContinuousParameter("x{0}".format(i), -1, 1) for i in range(1, 3)]) @@ -216,8 +224,8 @@ def test_constrained_EI(self): X = design.generate() Yo = parabola2d(X) Yc = plane(X) - m1 = GPflow.gpr.GPR(X, Yo, GPflow.kernels.RBF(2, ARD=True)) - m2 = GPflow.gpr.GPR(X, Yc, GPflow.kernels.RBF(2, ARD=True)) + m1 = GPflow.gpr.GPR(X, Yo, GPflow.kernels.RBF(2, ARD=True, lengthscales=X.std(axis=0))) + m2 = GPflow.gpr.GPR(X, Yc, GPflow.kernels.RBF(2, ARD=True, lengthscales=X.std(axis=0) / 2)) joint = GPflowOpt.acquisition.ExpectedImprovement(m1) * GPflowOpt.acquisition.ProbabilityOfFeasibility(m2) np.testing.assert_allclose(joint.objective_indices(), np.array([0], dtype=int)) @@ -228,10 +236,12 @@ def test_hierarchy(self): X = design.generate() Yo = parabola2d(X) Yc = plane(X) - m1 = GPflow.gpr.GPR(X, Yo, GPflow.kernels.RBF(2, ARD=True)) - m2 = GPflow.gpr.GPR(X, Yc, GPflow.kernels.RBF(2, ARD=True)) + m1 = GPflow.gpr.GPR(X, Yo, GPflow.kernels.RBF(2, ARD=True, lengthscales=X.std(axis=0))) + m2 = GPflow.gpr.GPR(X, Yo, GPflow.kernels.RBF(2, ARD=True, lengthscales=X.std(axis=0))) + m3 = GPflow.gpr.GPR(X, Yc, GPflow.kernels.RBF(2, ARD=True, lengthscales=X.std(axis=0) / 2)) joint = GPflowOpt.acquisition.ExpectedImprovement(m1) * \ - (GPflowOpt.acquisition.ProbabilityOfFeasibility(m2) + GPflowOpt.acquisition.ExpectedImprovement(m1)) + (GPflowOpt.acquisition.ProbabilityOfFeasibility(m3) + + GPflowOpt.acquisition.ExpectedImprovement(m2)) np.testing.assert_allclose(joint.objective_indices(), np.array([0, 2], dtype=int)) np.testing.assert_allclose(joint.constraint_indices(), np.array([1], dtype=int)) diff --git a/testing/test_optimizers.py b/testing/test_optimizers.py index 0b76eba..69420cb 100644 --- a/testing/test_optimizers.py +++ b/testing/test_optimizers.py @@ -134,7 +134,7 @@ def setUp(self): super(TestBayesianOptimizer, self).setUp() design = GPflowOpt.design.FactorialDesign(4, self.domain) X, Y = design.generate(), parabola2d(design.generate())[0] - model = GPflow.gpr.GPR(X, Y, GPflow.kernels.RBF(2, ARD=True)) + model = GPflow.gpr.GPR(X, Y, GPflow.kernels.RBF(2, ARD=True, lengthscales=X.std(axis=0))) acquisition = GPflowOpt.acquisition.ExpectedImprovement(model) self.optimizer = GPflowOpt.BayesianOptimizer(self.domain, acquisition) From 58a41d3ead7475b83e2a0c5551b0e8635607975d Mon Sep 17 00:00:00 2001 From: Joachim van der Herten Date: Fri, 5 May 2017 11:32:10 +0200 Subject: [PATCH 2/5] Fixed lengthscales for python 2.7 --- testing/test_optimizers.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/testing/test_optimizers.py b/testing/test_optimizers.py index 69420cb..347c2df 100644 --- a/testing/test_optimizers.py +++ b/testing/test_optimizers.py @@ -134,7 +134,7 @@ def setUp(self): super(TestBayesianOptimizer, self).setUp() design = GPflowOpt.design.FactorialDesign(4, self.domain) X, Y = design.generate(), parabola2d(design.generate())[0] - model = GPflow.gpr.GPR(X, Y, GPflow.kernels.RBF(2, ARD=True, lengthscales=X.std(axis=0))) + model = GPflow.gpr.GPR(X, Y, GPflow.kernels.RBF(2, ARD=True, lengthscales=X.std(axis=0) / 4)) acquisition = GPflowOpt.acquisition.ExpectedImprovement(model) self.optimizer = GPflowOpt.BayesianOptimizer(self.domain, acquisition) From 95d5abe7d12b7952ba05a990c91f13ff71bfc038 Mon Sep 17 00:00:00 2001 From: Joachim van der Herten Date: Fri, 5 May 2017 12:22:44 +0200 Subject: [PATCH 3/5] Prior should solve the optimization instabilities --- testing/test_acquisition.py | 4 +++- testing/test_optimizers.py | 3 ++- 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/testing/test_acquisition.py b/testing/test_acquisition.py index 84af083..d5a692c 100644 --- a/testing/test_acquisition.py +++ b/testing/test_acquisition.py @@ -26,7 +26,9 @@ def create_parabola_model(self, design=None): if design is None: design = GPflowOpt.design.FactorialDesign(4, self.domain) X, Y = design.generate(), parabola2d(design.generate()) - m = GPflow.gpr.GPR(X, Y, GPflow.kernels.RBF(2, ARD=True, lengthscales=X.std(axis=0))) + m = GPflow.gpr.GPR(X, Y, GPflow.kernels.RBF(2, ARD=True, lengthscales=X.var(axis=0))) + m.kern.variance.prior = GPflow.priors.Gamma(3,1/3) + print(m) return m def create_plane_model(self, design=None): diff --git a/testing/test_optimizers.py b/testing/test_optimizers.py index 347c2df..9375d44 100644 --- a/testing/test_optimizers.py +++ b/testing/test_optimizers.py @@ -134,7 +134,8 @@ def setUp(self): super(TestBayesianOptimizer, self).setUp() design = GPflowOpt.design.FactorialDesign(4, self.domain) X, Y = design.generate(), parabola2d(design.generate())[0] - model = GPflow.gpr.GPR(X, Y, GPflow.kernels.RBF(2, ARD=True, lengthscales=X.std(axis=0) / 4)) + model = GPflow.gpr.GPR(X, Y, GPflow.kernels.RBF(2, ARD=True, lengthscales=X.var(axis=0))) + model.kern.variance.prior = GPflow.priors.Gamma(3, 1 / 3) acquisition = GPflowOpt.acquisition.ExpectedImprovement(model) self.optimizer = GPflowOpt.BayesianOptimizer(self.domain, acquisition) From dd5aea49c9434c4d59c3b9df99bc3dd406426b3d Mon Sep 17 00:00:00 2001 From: Joachim van der Herten Date: Fri, 5 May 2017 17:01:46 +0200 Subject: [PATCH 4/5] Implemented a restart strategy for optimization of the hypers --- GPflowOpt/acquisition.py | 25 +++++++++++++++++-------- testing/test_acquisition.py | 5 ++--- 2 files changed, 19 insertions(+), 11 deletions(-) diff --git a/GPflowOpt/acquisition.py b/GPflowOpt/acquisition.py index e8f6c21..adcd209 100644 --- a/GPflowOpt/acquisition.py +++ b/GPflowOpt/acquisition.py @@ -37,21 +37,30 @@ class Acquisition(Parameterized): optimization) """ - def __init__(self, models=[]): + def __init__(self, models=[], optimize_restarts=5): super(Acquisition, self).__init__() self.models = ParamList(np.atleast_1d(models).tolist()) self._default_params = map(lambda m: m.get_free_state(), self.models) + + assert(optimize_restarts >= 0) + self._optimize_restarts = optimize_restarts self._optimize_all() def _optimize_all(self): for model, hypers in zip(self.models, self._default_params): - try: - model.optimize() - except: - # After data update, the starting point for the hypers may result in non-invertible K - # Reset the hypers to the values when the acquisition function was initialized. - model.set_state(hypers) - model.optimize() + runs = [] + # Start from supplied hyperparameters + model.set_state(hypers) + for i in range(self._optimize_restarts): + if i > 0: + model.randomize() + try: + result = model.optimize() + runs.append(result) + except tf.errors.InvalidArgumentError: + print("Warning - optimization restart {0}/{1} failed".format(i + 1, self._optimize_restarts)) + best_idx = np.argmin(map(lambda r: r.fun, runs)) + model.set_state(runs[best_idx].x) def _build_acquisition_wrapper(self, Xcand, gradients=True): acq = self.build_acquisition(Xcand) diff --git a/testing/test_acquisition.py b/testing/test_acquisition.py index d5a692c..25a5e2d 100644 --- a/testing/test_acquisition.py +++ b/testing/test_acquisition.py @@ -27,8 +27,7 @@ def create_parabola_model(self, design=None): design = GPflowOpt.design.FactorialDesign(4, self.domain) X, Y = design.generate(), parabola2d(design.generate()) m = GPflow.gpr.GPR(X, Y, GPflow.kernels.RBF(2, ARD=True, lengthscales=X.var(axis=0))) - m.kern.variance.prior = GPflow.priors.Gamma(3,1/3) - print(m) + m.kern.variance.prior = GPflow.priors.Gamma(3.0,1.0/3.0) return m def create_plane_model(self, design=None): @@ -107,7 +106,7 @@ def test_setup(self): p = np.array([[0.0, 0.0]]) self.acquisition.set_data(np.vstack((self.model.X.value, p)), np.vstack((self.model.Y.value, parabola2d(p)))) - self.assertTrue(np.allclose(self.acquisition.fmin.value, 0, atol=1e-2), msg="fmin not updated") + self.assertTrue(np.allclose(self.acquisition.fmin.value, 0, atol=1e-1), msg="fmin not updated") def test_EI_validity(self): Xcenter = np.random.rand(20, 2) * 0.25 - 0.125 From d69c046421768ec3af7a8062d4090990c7c6d4c1 Mon Sep 17 00:00:00 2001 From: Joachim van der Herten Date: Fri, 5 May 2017 17:14:40 +0200 Subject: [PATCH 5/5] Fixed initialization error of Gamma prior for python 2.7 --- testing/test_acquisition.py | 2 +- testing/test_optimizers.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/testing/test_acquisition.py b/testing/test_acquisition.py index 25a5e2d..e6afb6d 100644 --- a/testing/test_acquisition.py +++ b/testing/test_acquisition.py @@ -27,7 +27,7 @@ def create_parabola_model(self, design=None): design = GPflowOpt.design.FactorialDesign(4, self.domain) X, Y = design.generate(), parabola2d(design.generate()) m = GPflow.gpr.GPR(X, Y, GPflow.kernels.RBF(2, ARD=True, lengthscales=X.var(axis=0))) - m.kern.variance.prior = GPflow.priors.Gamma(3.0,1.0/3.0) + m.kern.variance.prior = GPflow.priors.Gamma(3.0, 1.0 / 3.0) return m def create_plane_model(self, design=None): diff --git a/testing/test_optimizers.py b/testing/test_optimizers.py index 9375d44..5f6cdf5 100644 --- a/testing/test_optimizers.py +++ b/testing/test_optimizers.py @@ -135,7 +135,7 @@ def setUp(self): design = GPflowOpt.design.FactorialDesign(4, self.domain) X, Y = design.generate(), parabola2d(design.generate())[0] model = GPflow.gpr.GPR(X, Y, GPflow.kernels.RBF(2, ARD=True, lengthscales=X.var(axis=0))) - model.kern.variance.prior = GPflow.priors.Gamma(3, 1 / 3) + model.kern.variance.prior = GPflow.priors.Gamma(3, 1.0 / 3.0) acquisition = GPflowOpt.acquisition.ExpectedImprovement(model) self.optimizer = GPflowOpt.BayesianOptimizer(self.domain, acquisition)