Handle deterministic predictions in Expected Improvement

bayes_opt/acquisition.py

@@ -845,8 +845,10 @@ def base_acq(self, mean: NDArray[Float], std: NDArray[Float]) -> NDArray[Float]:
             )
             raise ValueError(msg)
         a = mean - self.y_max - self.xi
-        z = a / std
-        return a * norm.cdf(z) + std * norm.pdf(z)
+        z = a / np.where(std == 0, 1, std)
+        improvement = a * norm.cdf(z) + std * norm.pdf(z)
+        # A deterministic prediction has no uncertainty to integrate over.
+        return np.where(std == 0, np.maximum(a, 0), improvement)
 
     def suggest(
         self,

tests/test_acquisition.py

@@ -7,7 +7,9 @@
 import pytest
 from scipy.optimize import NonlinearConstraint
 from scipy.spatial.distance import pdist
+from scipy.stats import norm
 from sklearn.gaussian_process import GaussianProcessRegressor
+from sklearn.gaussian_process.kernels import ConstantKernel
 
 from bayes_opt import BayesianOptimization, acquisition, exception
 from bayes_opt.acquisition import (
@@ -238,6 +240,27 @@ def test_expected_improvement(gp, target_space, random_state):
     assert acq.xi == 0.01
 
 
+@pytest.mark.parametrize("xi", [0.0, 0.25])
+def test_expected_improvement_zero_variance(xi):
+    acq = ExpectedImprovement(xi=xi)
+    acq.y_max = 1.0
+    mean = np.array([0.0, 1.0 + xi, 2.0, 1.5])
+    expected = [0.0, 0.0, 1.0 - xi, norm.expect(lambda x: max(x - 1.0 - xi, 0), loc=1.5, scale=0.5)]
+    with np.errstate(divide="raise", invalid="raise"):
+        np.testing.assert_allclose(acq.base_acq(mean, np.array([0.0, 0.0, 0.0, 0.5])), expected)
+        assert acq.base_acq(1.0 + xi, 0.0) == 0.0
+
+
+def test_expected_improvement_deterministic_gp():
+    gp = GaussianProcessRegressor(kernel=ConstantKernel(1.0, "fixed"), alpha=0, optimizer=None)
+    gp.fit([[0.0]], [1.0])
+    mean, std = gp.predict([[0.0], [1.0]], return_std=True)
+    np.testing.assert_array_equal(std, 0.0)
+    acq = ExpectedImprovement(xi=0.0)
+    acq.y_max = 1.0
+    np.testing.assert_array_equal(acq.base_acq(mean, std), 0.0)
+
+
 def test_expected_improvement_with_constraints(gp, constrained_target_space, random_state):
     acq = acquisition.ExpectedImprovement(exploration_decay=0.5, exploration_decay_delay=2, xi=0.01)
     assert acq.xi == 0.01