Repository: bayesian-optimization/BayesianOptimization

Issue:
Next point suggestion with mixed-integer optimisation
Hi all,

I see that mixed-integer optimisation has been allowed in the library. This is a very good improvement. Looking at the code, I see that the suggestion of the next point still happens with L-BFGS-B, which is meant for continuous optimisation. Is that the case? Or is there a part of the code I missed?

If so, I propose to implement something able to an optimiser able to deal with mixed-integer. I propose applying the differential evolution with L-BFGS-B refinement on the continuous variables algorithm that is implemented in scipy. However, I am open to suggestions.

Is this a wanted feature for this project? If so, I can work on the implementation.

Discussion:
@till-m: Hi @udicaprio,

good catch! However, the L-BFGSB optimization is indeed only over continuous parameters, see the relevant lines starting [roughly here](https://github.com/bayesian-optimization/BayesianOptimization/blob/232f9d95ce17aeeb6bfddb884d62c65e41292c98/bayes_opt/acquisition.py#L298).
In short, it does L-BFGS-B optimization over the continuous dimensions, starting from the `n` best randomly sampled points. This might not be the ideal way of doing it, in the paper they don't really discuss the optimization of the acquisition function very much. In my (very limited) experiments, this worked reasonably well, but I did feel it could be a bit finicky. Your approach sounds like it could be significantly more stable.

In this sense: I think this could really improve the library!@till-m: If you do decide to work on this, I think it would make sense to hold off on the full 3.0 release for now, depending on how long you think it will take you (no pressure though).@udicaprio: Hi @till-m,

Thank you. I will be working on this in the coming weeks. 
What is the deadline for the full release of 3.0?@till-m: Hey @udicaprio,

nothing fixed, really. A few weeks is great honestly! @udicaprio: The pull request implementing these functions has been merged

Solution patch:
diff --git a/bayes_opt/acquisition.py b/bayes_opt/acquisition.py
index a81a6df..c29b117 100644
--- a/bayes_opt/acquisition.py
+++ b/bayes_opt/acquisition.py
@@ -27,7 +27,9 @@ from typing import TYPE_CHECKING, Any, Literal, NoReturn
 
 import numpy as np
 from numpy.random import RandomState
-from scipy.optimize import minimize
+from packaging import version
+from scipy import __version__ as scipy_version
+from scipy.optimize._differentialevolution import DifferentialEvolutionSolver, minimize
 from scipy.special import softmax
 from scipy.stats import norm
 from sklearn.gaussian_process import GaussianProcessRegressor
@@ -38,6 +40,7 @@ from bayes_opt.exception import (
     TargetSpaceEmptyError,
 )
 from bayes_opt.target_space import TargetSpace
+from bayes_opt.util import ensure_rng
 
 if TYPE_CHECKING:
     from collections.abc import Callable, Sequence
@@ -60,13 +63,7 @@ class AcquisitionFunction(abc.ABC):
     """
 
     def __init__(self, random_state: int | RandomState | None = None) -> None:
-        if random_state is not None:
-            if isinstance(random_state, RandomState):
-                self.random_state = random_state
-            else:
-                self.random_state = RandomState(random_state)
-        else:
-            self.random_state = RandomState()
+        self.random_state = ensure_rng(random_state)
         self.i = 0
 
     def _serialize_random_state(self) -> dict | None:
@@ -274,7 +271,7 @@ class AcquisitionFunction(abc.ABC):
             acq, space, n_random=max(n_random, n_l_bfgs_b), n_x_seeds=n_l_bfgs_b
         )
         if n_l_bfgs_b:
-            x_min_l, min_acq_l = self._l_bfgs_b_minimize(acq, space, x_seeds=x_seeds)
+            x_min_l, min_acq_l = self._smart_minimize(acq, space, x_seeds=x_seeds)
             # Either n_random or n_l_bfgs_b is not 0 => at least one of x_min_r and x_min_l is not None
             if min_acq_r > min_acq_l:
                 return x_min_l
@@ -323,7 +320,7 @@ class AcquisitionFunction(abc.ABC):
             x_seeds = []
         return x_min, min_acq, x_seeds
 
-    def _l_bfgs_b_minimize(
+    def _smart_minimize(
         self,
         acq: Callable[[NDArray[Float]], NDArray[Float]],
         space: TargetSpace,
@@ -353,38 +350,66 @@ class AcquisitionFunction(abc.ABC):
         continuous_dimensions = space.continuous_dimensions
         continuous_bounds = space.bounds[continuous_dimensions]
 
-        if not continuous_dimensions.any():
-            min_acq = np.inf
-            x_min = np.array([np.nan] * space.bounds.shape[0])
-            return x_min, min_acq
-
         min_acq: float | None = None
         x_try: NDArray[Float]
         x_min: NDArray[Float]
-        for x_try in x_seeds:
 
-            def continuous_acq(x: NDArray[Float], x_try=x_try) -> NDArray[Float]:
-                x_try[continuous_dimensions] = x
-                return acq(x_try)
+        # Case of continous optimization
+        if all(continuous_dimensions):
+            for x_try in x_seeds:
+                res: OptimizeResult = minimize(acq, x_try, bounds=continuous_bounds, method="L-BFGS-B")
+                if not res.success:
+                    continue
 
-            # Find the minimum of minus the acquisition function
-            res: OptimizeResult = minimize(
-                continuous_acq, x_try[continuous_dimensions], bounds=continuous_bounds, method="L-BFGS-B"
-            )
-            # See if success
-            if not res.success:
-                continue
+                # Store it if better than previous minimum(maximum).
+                if min_acq is None or np.squeeze(res.fun) >= min_acq:
+                    x_try = res.x
+                    x_min = x_try
+                    min_acq = np.squeeze(res.fun)
 
-            # Store it if better than previous minimum(maximum).
-            if min_acq is None or np.squeeze(res.fun) >= min_acq:
-                x_try[continuous_dimensions] = res.x
-                x_min = x_try
-                min_acq = np.squeeze(res.fun)
+        # Case of mixed-integer optimization
+        else:
+            ntrials = max(1, len(x_seeds) // 100)
+
+            for _ in range(ntrials):
+                xinit = space.random_sample(15 * len(space.bounds), random_state=self.random_state)
+
+                de_parameters = {"func": acq, "bounds": space.bounds, "polish": False, "init": xinit}
+                if version.parse(scipy_version) < version.parse("1.15.0"):
+                    de_parameters["seed"] = self.random_state
+                else:
+                    de_parameters["rng"] = self.random_state
+
+                de = DifferentialEvolutionSolver(**de_parameters)
+                res_de: OptimizeResult = de.solve()
+                # Check if success
+                if not res_de.success:
+                    continue
+
+                x_min = res_de.x
+                min_acq = np.squeeze(res_de.fun)
+
+                # Refine the identification of continous parameters with deterministic search
+                if any(continuous_dimensions):
+                    x_try = x_min.copy()
+
+                    def continuous_acq(x: NDArray[Float], x_try=x_try) -> NDArray[Float]:
+                        x_try[continuous_dimensions] = x
+                        return acq(x_try)
+
+                    res: OptimizeResult = minimize(
+                        continuous_acq, x_min[continuous_dimensions], bounds=continuous_bounds
+                    )
+                    if np.squeeze(res.fun) >= min_acq and res.success:
+                        x_try[continuous_dimensions] = res.x
+                        x_min = x_try
+                        min_acq = np.squeeze(res.fun)
 
         if min_acq is None:
             min_acq = np.inf
             x_min = np.array([np.nan] * space.bounds.shape[0])
 
+        x_min = space.kernel_transform(x_min).reshape(x_min.shape)
         # Clip output to make sure it lies within the bounds. Due to floating
         # point technicalities this is not always the case.
         return np.clip(x_min, space.bounds[:, 0], space.bounds[:, 1]), min_acq
diff --git a/bayes_opt/target_space.py b/bayes_opt/target_space.py
index 39d1f99..ba41bbf 100644
--- a/bayes_opt/target_space.py
+++ b/bayes_opt/target_space.py
@@ -312,9 +312,7 @@ class TargetSpace:
             Representation of the parameters as an array.
         """
         if set(params) != set(self.keys):
-            error_msg = (
-                f"Parameters' keys ({params}) do " f"not match the expected set of keys ({self.keys})."
-            )
+            error_msg = f"Parameters' keys ({params}) do not match the expected set of keys ({self.keys})."
             raise ValueError(error_msg)
         return self._to_float(params)
 
@@ -356,15 +354,14 @@ class TargetSpace:
         """
         if len(x) != self._dim:
             error_msg = (
-                f"Size of array ({len(x)}) is different than the "
-                f"expected number of parameters ({self._dim})."
+                f"Size of array ({len(x)}) is different than the expected number of parameters ({self._dim})."
             )
             raise ValueError(error_msg)
         return self._to_params(x)
 
     def _to_float(self, value: Mapping[str, float | NDArray[Float]]) -> NDArray[Float]:
         if set(value) != set(self.keys):
-            msg = f"Parameters' keys ({value}) do " f"not match the expected set of keys ({self.keys})."
+            msg = f"Parameters' keys ({value}) do not match the expected set of keys ({self.keys})."
             raise ValueError(msg)
         res = np.zeros(self._dim)
         for key in self._keys:
@@ -704,9 +701,7 @@ class TargetSpace:
                 params_config[key] = new_params_config[key]
             dims = dims + params_config[key].dim
         if dims != self.dim:
-            msg = (
-                f"Dimensions of new bounds ({dims}) does not match" f" dimensions of old bounds ({self.dim})."
-            )
+            msg = f"Dimensions of new bounds ({dims}) does not match dimensions of old bounds ({self.dim})."
             raise ValueError(msg)
         self._params_config = params_config
         self._bounds = self.calculate_bounds()


Test patch:
diff --git a/tests/test_acquisition.py b/tests/test_acquisition.py
index e3d3e51..0e46aba 100644
--- a/tests/test_acquisition.py
+++ b/tests/test_acquisition.py
@@ -140,7 +140,7 @@ def test_upper_confidence_bound(gp, target_space, random_state):
     assert acq.kappa == 0.5
 
 
-def test_l_bfgs_fails(target_space, random_state):
+def test_smart_minimize_fails(target_space, random_state):
     acq = acquisition.UpperConfidenceBound(random_state=random_state)
 
     def fun(x):
@@ -149,7 +149,7 @@ def test_l_bfgs_fails(target_space, random_state):
         except IndexError:
             return np.nan
 
-    _, min_acq_l = acq._l_bfgs_b_minimize(fun, space=target_space, x_seeds=np.array([[2.5, 0.5]]))
+    _, min_acq_l = acq._smart_minimize(fun, space=target_space, x_seeds=np.array([[2.5, 0.5]]))
     assert min_acq_l == np.inf