diff --git a/AGENTS.md b/AGENTS.md
index 59629ab1fb..0696db1519 100644
--- a/AGENTS.md
+++ b/AGENTS.md
@@ -105,7 +105,7 @@ Custom `@classproperty` from `baybe.utils.basic` for class-level computed proper
 | Element | Convention | Examples |
 |---------|------------|---------|
 | Variables/functions | `snake_case` | `batch_size`, `add_measurements` |
-| Classes | `PascalCase` | `Campaign`, `BotorchRecommender` |
+| Classes | `PascalCase` | `Campaign`, `BayesianRecommender` |
 | Constants | `SCREAMING_SNAKE_CASE` | `_RECOMMENDED`, `_TYPE_FIELD` |
 | Private members | `_` prefix | `_cached_recommendation`, `_fit()` |
 | Booleans | `is_`/`has_`/`supports_` | `is_numerical`, `supports_transfer_learning` |
diff --git a/README.md b/README.md
index 6cfe839f20..064a72bc75 100644
--- a/README.md
+++ b/README.md
@@ -234,14 +234,14 @@ For our example, we combine two recommenders via a so-called meta recommender na
 
 ```python
 from baybe.recommenders import (
-    BotorchRecommender,
+    BayesianRecommender,
     FPSRecommender,
     TwoPhaseMetaRecommender,
 )
 
 recommender = TwoPhaseMetaRecommender(
     initial_recommender=FPSRecommender(),  # farthest point sampling
-    recommender=BotorchRecommender(),  # Bayesian model-based optimization
+    recommender=BayesianRecommender(),  # Bayesian model-based optimization
 )
 ```
 
diff --git a/baybe/recommenders/__init__.py b/baybe/recommenders/__init__.py
index 88bc984337..e35da9b05b 100644
--- a/baybe/recommenders/__init__.py
+++ b/baybe/recommenders/__init__.py
@@ -6,6 +6,7 @@
     TwoPhaseMetaRecommender,
 )
 from baybe.recommenders.naive import NaiveHybridSpaceRecommender
+from baybe.recommenders.pure.bayesian.base import BayesianRecommender
 from baybe.recommenders.pure.bayesian.botorch import BotorchRecommender
 from baybe.recommenders.pure.nonpredictive.clustering import (
     GaussianMixtureClusteringRecommender,
@@ -18,6 +19,7 @@
 )
 
 __all__ = [
+    "BayesianRecommender",
     "BotorchRecommender",
     "FPSRecommender",
     "GaussianMixtureClusteringRecommender",
diff --git a/baybe/recommenders/meta/sequential.py b/baybe/recommenders/meta/sequential.py
index 8f581c74e0..1c560d6f35 100644
--- a/baybe/recommenders/meta/sequential.py
+++ b/baybe/recommenders/meta/sequential.py
@@ -18,7 +18,7 @@
 from baybe.objectives.base import Objective
 from baybe.recommenders.base import RecommenderProtocol
 from baybe.recommenders.meta.base import MetaRecommender
-from baybe.recommenders.pure.bayesian.botorch import BotorchRecommender
+from baybe.recommenders.pure.bayesian import BayesianRecommender
 from baybe.recommenders.pure.nonpredictive.sampling import RandomRecommender
 from baybe.searchspace import SearchSpace
 from baybe.serialization import (
@@ -52,7 +52,7 @@ class TwoPhaseMetaRecommender(MetaRecommender):
     """The initial recommender used by the meta recommender."""
 
     recommender: RecommenderProtocol = field(
-        factory=BotorchRecommender, validator=instance_of(RecommenderProtocol)
+        factory=BayesianRecommender, validator=instance_of(RecommenderProtocol)
     )
     """The recommender used by the meta recommender after the switch."""
 
diff --git a/baybe/recommenders/naive.py b/baybe/recommenders/naive.py
index 5b602d881b..7afefcd76f 100644
--- a/baybe/recommenders/naive.py
+++ b/baybe/recommenders/naive.py
@@ -10,7 +10,6 @@
 from baybe.objectives.base import Objective
 from baybe.recommenders.pure.base import PureRecommender
 from baybe.recommenders.pure.bayesian.base import BayesianRecommender
-from baybe.recommenders.pure.bayesian.botorch import BotorchRecommender
 from baybe.recommenders.pure.nonpredictive.base import NonPredictiveRecommender
 from baybe.searchspace import SearchSpace, SearchSpaceType
 from baybe.utils.dataframe import to_tensor
@@ -39,13 +38,13 @@ class NaiveHybridSpaceRecommender(PureRecommender):
     # works for now. Still, we manually check whether the disc_recommender belongs to
     # one of these two subclasses such that we might be able to easily spot a potential
     # problem that might come up when implementing new subclasses of PureRecommender
-    disc_recommender: PureRecommender = field(factory=BotorchRecommender)
+    disc_recommender: PureRecommender = field(factory=BayesianRecommender)
     """The recommender used for the discrete subspace. Default:
-    :class:`baybe.recommenders.pure.bayesian.botorch.core.BotorchRecommender`"""
+    :class:`baybe.recommenders.pure.bayesian.base.BayesianRecommender`"""
 
-    cont_recommender: BayesianRecommender = field(factory=BotorchRecommender)
+    cont_recommender: BayesianRecommender = field(factory=BayesianRecommender)
     """The recommender used for the continuous subspace. Default:
-    :class:`baybe.recommenders.pure.bayesian.botorch.core.BotorchRecommender`"""
+    :class:`baybe.recommenders.pure.bayesian.base.BayesianRecommender`"""
 
     @override
     def recommend(
diff --git a/baybe/recommenders/pure/__init__.py b/baybe/recommenders/pure/__init__.py
index 76340a57ad..1b6d0581d4 100644
--- a/baybe/recommenders/pure/__init__.py
+++ b/baybe/recommenders/pure/__init__.py
@@ -4,6 +4,7 @@
 recommendations. They can be part of meta recommenders.
 """
 
+from baybe.recommenders.pure.bayesian.base import BayesianRecommender
 from baybe.recommenders.pure.bayesian.botorch import BotorchRecommender
 from baybe.recommenders.pure.nonpredictive import (
     FPSRecommender,
@@ -14,6 +15,7 @@
 )
 
 __all__ = [
+    "BayesianRecommender",
     "BotorchRecommender",
     "FPSRecommender",
     "GaussianMixtureClusteringRecommender",
diff --git a/baybe/recommenders/pure/bayesian/__init__.py b/baybe/recommenders/pure/bayesian/__init__.py
index 9fa0daf758..771cbf0fab 100644
--- a/baybe/recommenders/pure/bayesian/__init__.py
+++ b/baybe/recommenders/pure/bayesian/__init__.py
@@ -1,7 +1,9 @@
 """Bayesian recommenders."""
 
+from baybe.recommenders.pure.bayesian.base import BayesianRecommender
 from baybe.recommenders.pure.bayesian.botorch import BotorchRecommender
 
 __all__ = [
+    "BayesianRecommender",
     "BotorchRecommender",
 ]
diff --git a/baybe/recommenders/pure/bayesian/base.py b/baybe/recommenders/pure/bayesian/base.py
index db163a8556..8273fbb454 100644
--- a/baybe/recommenders/pure/bayesian/base.py
+++ b/baybe/recommenders/pure/bayesian/base.py
@@ -4,11 +4,14 @@
 
 import gc
 from abc import ABC
-from typing import TYPE_CHECKING
+from collections.abc import Callable, Iterable
+from typing import TYPE_CHECKING, Any
 
+import numpy as np
 import pandas as pd
 from attrs import define, field
 from attrs.converters import optional
+from attrs.validators import ge, gt, instance_of
 from typing_extensions import override
 
 from baybe.acquisition import qLogEI, qLogNEHVI
@@ -16,20 +19,40 @@
 from baybe.acquisition.utils import convert_acqf
 from baybe.exceptions import (
     IncompatibleAcquisitionFunctionError,
+    InfeasibilityError,
 )
 from baybe.objectives.base import Objective
 from baybe.recommenders.pure.base import PureRecommender
-from baybe.searchspace import SearchSpace
+from baybe.recommenders.pure.bayesian.botorch.optimizers.base import OptimizerProtocol
+from baybe.recommenders.pure.bayesian.botorch.optimizers.basic import GradientOptimizer
+from baybe.recommenders.pure.bayesian.continuous import (
+    recommend_continuous_torch,
+)
+from baybe.recommenders.pure.bayesian.discrete import (
+    recommend_discrete_with_subsets,
+    recommend_discrete_without_subsets,
+)
+from baybe.recommenders.pure.bayesian.hybrid import (
+    recommend_hybrid_with_subsets,
+    recommend_hybrid_without_subsets,
+)
+from baybe.searchspace import (
+    SearchSpace,
+    SubspaceContinuous,
+    SubspaceDiscrete,
+)
 from baybe.settings import Settings
 from baybe.surrogates import GaussianProcessSurrogate
 from baybe.surrogates.base import (
     Surrogate,
     SurrogateProtocol,
 )
+from baybe.utils.sampling_algorithms import DiscreteSamplingMethod
 from baybe.utils.validation import preprocess_dataframe, validate_object_names
 
 if TYPE_CHECKING:
     from botorch.acquisition import AcquisitionFunction as BoAcquisitionFunction
+    from torch import Tensor
 
 
 def _autoreplicate(surrogate: SurrogateProtocol, /) -> SurrogateProtocol:
@@ -55,6 +78,39 @@ class BayesianRecommender(PureRecommender, ABC):
     )
     """The acquisition function. When omitted, a default is used."""
 
+    optimizer: OptimizerProtocol = field(
+        alias="optimizer",
+        default=GradientOptimizer(),
+    )
+    """The acquisition function optimizer."""
+
+    # TODO: Move fields to respective optimizers
+    hybrid_sampler: DiscreteSamplingMethod | None = field(
+        converter=optional(DiscreteSamplingMethod), default=None
+    )
+    """Strategy used for sampling the discrete subspace when performing hybrid search
+    space optimization."""
+
+    sampling_percentage: float = field(default=1.0)
+    """Percentage of discrete search space that is sampled when performing hybrid search
+    space optimization. Ignored when ``hybrid_sampler="None"``."""
+
+    n_restarts: int = field(validator=[instance_of(int), gt(0)], default=10)
+    """Number of times gradient-based optimization is restarted from different initial
+    points. **Does not affect purely discrete optimization**.
+    """
+
+    n_raw_samples: int = field(validator=[instance_of(int), gt(0)], default=64)
+    """Number of raw samples drawn for the initialization heuristic in gradient-based
+    optimization. **Does not affect purely discrete optimization**.
+    """
+
+    max_n_subsets: int = field(default=10, validator=[instance_of(int), ge(1)])
+    """Maximum number of subsets to evaluate when subset-generating constraints are
+    present (e.g., continuous cardinality constraints). If the total number of
+    subsets exceeds this limit, a random subset of that size is sampled for
+    optimization instead of performing an exhaustive search."""
+
     # TODO: The objective is currently only required for validating the recommendation
     #   context. Once multi-target support is complete, we might want to refactor
     #   the validation mechanism, e.g. by
@@ -67,6 +123,20 @@ class BayesianRecommender(PureRecommender, ABC):
     _botorch_acqf = field(default=None, init=False, eq=False)
     """The induced BoTorch acquisition function."""
 
+    @sampling_percentage.validator
+    def _validate_percentage(  # noqa: DOC101, DOC103
+        self, _: Any, value: float
+    ) -> None:
+        """Validate that the given value is in fact a percentage.
+
+        Raises:
+            ValueError: If ``value`` is not between 0 and 1.
+        """
+        if not 0 <= value <= 1:
+            raise ValueError(
+                f"Hybrid sampling percentage needs to be between 0 and 1 but is {value}"
+            )
+
     def _get_acquisition_function(self, objective: Objective) -> AcquisitionFunction:
         """Select the appropriate default acquisition function for the given context."""
         if self.acquisition_function is None:
@@ -196,6 +266,147 @@ def recommend(
             else:
                 raise
 
+    @override
+    def _recommend_discrete(
+        self,
+        subspace_discrete: SubspaceDiscrete,
+        candidates_exp: pd.DataFrame,
+        batch_size: int,
+    ) -> pd.Index:
+        """Generate recommendations from a discrete search space.
+
+        Dispatches to the appropriate optimization routine depending on whether
+        subset constraints are present.
+
+        Args:
+            subspace_discrete: The discrete subspace from which to generate
+                recommendations.
+            candidates_exp: The experimental representation of all discrete candidate
+                points to be considered.
+            batch_size: The size of the recommendation batch.
+
+        Returns:
+            The dataframe indices of the recommended points in the provided
+            experimental representation.
+        """
+        if subspace_discrete.n_subsets > 0:
+            return recommend_discrete_with_subsets(
+                self, subspace_discrete, candidates_exp, batch_size
+            )
+        return recommend_discrete_without_subsets(
+            self, subspace_discrete, candidates_exp, batch_size
+        )
+
+    @override
+    def _recommend_continuous(
+        self,
+        subspace_continuous: SubspaceContinuous,
+        batch_size: int,
+    ) -> pd.DataFrame:
+        """Generate recommendations from a continuous search space.
+
+        Args:
+            subspace_continuous: The continuous subspace from which to generate
+                recommendations.
+            batch_size: The size of the recommendation batch.
+
+        Raises:
+            IncompatibleAcquisitionFunctionError: If a non-Monte Carlo acquisition
+                function is used with a batch size > 1.
+
+        Returns:
+            A dataframe containing the recommendations as individual rows.
+        """
+        assert self._objective is not None
+        if (
+            batch_size > 1
+            and not self._get_acquisition_function(self._objective).supports_batching
+        ):
+            raise IncompatibleAcquisitionFunctionError(
+                f"The '{self.__class__.__name__}' only works with Monte Carlo "
+                f"acquisition functions for batch sizes > 1."
+            )
+
+        points, _ = recommend_continuous_torch(self, subspace_continuous, batch_size)
+
+        return pd.DataFrame(points, columns=subspace_continuous.parameter_names)
+
+    @override
+    def _recommend_hybrid(
+        self,
+        searchspace: SearchSpace,
+        candidates_exp: pd.DataFrame,
+        batch_size: int,
+    ) -> pd.DataFrame:
+        """Generate recommendations from a hybrid search space.
+
+        Dispatches to the appropriate optimization routine depending on whether
+        subset constraints are present.
+
+        Args:
+            searchspace: The search space in which the recommendations should be made.
+            candidates_exp: The experimental representation of the candidates
+                of the discrete subspace.
+            batch_size: The size of the calculated batch.
+
+        Returns:
+            The recommended points.
+        """
+        if searchspace.n_subsets > 0:
+            return recommend_hybrid_with_subsets(
+                self, searchspace, candidates_exp, batch_size
+            )
+        return recommend_hybrid_without_subsets(
+            self, searchspace, candidates_exp, batch_size
+        )
+
+    def _optimize_over_subsets(
+        self,
+        subset_callables: Iterable[Callable[[], tuple[Any, Tensor]]],
+    ) -> tuple[Any, Tensor]:
+        """Optimize across subsets and return the result with the best acqf value.
+
+        Each callable performs optimization for one subset configuration and returns
+        a ``(result, acquisition_value)`` tuple. Subsets that raise
+        ``InfeasibilityError`` are silently skipped.
+
+        Args:
+            subset_callables: An iterable of zero-argument callables. Each callable
+                runs the optimization for one subset and returns
+                ``(result, acqf_value)``. It may raise ``InfeasibilityError`` if the
+                subset is infeasible.
+
+        Raises:
+            InfeasibilityError: If none of the subsets has a feasible solution.
+
+        Returns:
+            The result and acquisition value of the best subset.
+        """
+        from botorch.exceptions.errors import InfeasibilityError as BoInfeasibilityError
+
+        results_all: list = []
+        acqf_values_all: list[Tensor] = []
+
+        for optimize_fn in subset_callables:
+            try:
+                result, acqf_value = optimize_fn()
+                results_all.append(result)
+                acqf_values_all.append(acqf_value)
+            except (BoInfeasibilityError, InfeasibilityError):
+                pass
+
+        if not results_all:
+            raise InfeasibilityError(
+                "No feasible solution could be found. Potentially the specified "
+                "constraints are too restrictive, i.e. there may be too many "
+                "constraints or thresholds may have been set too tightly. "
+                "Consider relaxing the constraints to improve the chances "
+                "of finding a feasible solution."
+            )
+
+        best_idx = np.argmax(acqf_values_all)
+        return results_all[best_idx], acqf_values_all[best_idx]
+
     def acquisition_values(
         self,
         candidates: pd.DataFrame,
diff --git a/baybe/recommenders/pure/bayesian/botorch/core.py b/baybe/recommenders/pure/bayesian/botorch/core.py
index 7953d5ca74..205e8cb8c2 100644
--- a/baybe/recommenders/pure/bayesian/botorch/core.py
+++ b/baybe/recommenders/pure/bayesian/botorch/core.py
@@ -4,295 +4,95 @@
 
 import gc
 import warnings
-from collections.abc import Callable, Iterable
-from typing import TYPE_CHECKING, Any, ClassVar
+from typing import TYPE_CHECKING
 
-import numpy as np
-import pandas as pd
-from attrs import define, field
-from attrs.converters import optional as optional_c
-from attrs.validators import ge, gt, instance_of
-from typing_extensions import override
-
-from baybe.exceptions import (
-    IncompatibleAcquisitionFunctionError,
-    InfeasibilityError,
-)
-from baybe.recommenders.pure.bayesian.base import BayesianRecommender
-from baybe.recommenders.pure.bayesian.botorch.continuous import (
-    recommend_continuous_torch,
-)
-from baybe.recommenders.pure.bayesian.botorch.discrete import (
-    recommend_discrete_with_subsets,
-    recommend_discrete_without_subsets,
-)
-from baybe.recommenders.pure.bayesian.botorch.hybrid import (
-    recommend_hybrid_with_subsets,
-    recommend_hybrid_without_subsets,
-)
-from baybe.searchspace import (
-    SearchSpace,
-    SearchSpaceType,
-    SubspaceContinuous,
-    SubspaceDiscrete,
-)
-from baybe.utils.conversion import to_string
+from baybe.acquisition.base import AcquisitionFunction
+from baybe.recommenders.pure.bayesian.botorch.optimizers.basic import GradientOptimizer
+from baybe.surrogates.base import SurrogateProtocol
+from baybe.surrogates.gaussian_process.core import GaussianProcessSurrogate
 from baybe.utils.sampling_algorithms import DiscreteSamplingMethod
 
 if TYPE_CHECKING:
-    from torch import Tensor
-
-
-@define(kw_only=True)
-class BotorchRecommender(BayesianRecommender):
-    """A pure recommender utilizing Botorch's optimization machinery.
-
-    This recommender makes use of Botorch's ``optimize_acqf_discrete``,
-    ``optimize_acqf`` and ``optimize_acqf_mixed`` functions to optimize discrete,
-    continuous and hybrid search spaces, respectively. Accordingly, it can be applied to
-    all kinds of search spaces.
-
-    Note:
-        In hybrid search spaces, the used algorithm performs a brute-force optimization
-        that can be computationally expensive. Thus, the behavior of the algorithm in
-        hybrid search spaces can be controlled via two additional parameters.
-    """
-
-    # Class variables
-    compatibility: ClassVar[SearchSpaceType] = SearchSpaceType.HYBRID
-    # See base class.
-
-    supports_discrete_subset_generating_constraints: ClassVar[bool] = True
-    # See base class.
-
-    # Object variables
-    sequential_continuous: bool = field(default=True)
-    """Flag defining whether to apply sequential greedy or batch optimization in
-    **continuous** search spaces. In discrete/hybrid spaces, sequential greedy
-    optimization is applied automatically.
-    """
-
-    hybrid_sampler: DiscreteSamplingMethod | None = field(
-        converter=optional_c(DiscreteSamplingMethod), default=None
-    )
-    """Strategy used for sampling the discrete subspace when performing hybrid search
-    space optimization."""
-
-    sampling_percentage: float = field(default=1.0)
-    """Percentage of discrete search space that is sampled when performing hybrid search
-    space optimization. Ignored when ``hybrid_sampler="None"``."""
-
-    n_restarts: int = field(validator=[instance_of(int), gt(0)], default=10)
-    """Number of times gradient-based optimization is restarted from different initial
-    points. **Does not affect purely discrete optimization**.
-    """
-
-    n_raw_samples: int = field(validator=[instance_of(int), gt(0)], default=64)
-    """Number of raw samples drawn for the initialization heuristic in gradient-based
-    optimization. **Does not affect purely discrete optimization**.
+    from baybe.recommenders.pure.bayesian.base import BayesianRecommender
+
+
+def BotorchRecommender(
+    *,
+    surrogate_model: SurrogateProtocol = GaussianProcessSurrogate(),
+    acquisition_function: AcquisitionFunction | None = None,
+    sequential_continuous: bool = True,
+    hybrid_sampler: DiscreteSamplingMethod | None = None,
+    sampling_percentage: float = 1.0,
+    n_restarts: int = 10,
+    n_raw_samples: int = 64,
+    max_n_subsets: int = 10,
+    max_n_subspaces: int | None = None,
+) -> BayesianRecommender:
+    """Use factory function for BotorchRecommender deprecation.
+
+    This recommender will be deprecated in a future version.
+    This function provides the interface for creating a BayesianRecommender
+    based on a BotorchRecommender.
+
+    Args:
+        surrogate_model: The surrogate model to be used.
+        acquisition_function: The acquisition function to be used.
+        sequential_continuous: See :class:`BayesianRecommender`.
+        hybrid_sampler: See :class:`BayesianRecommender`.
+        sampling_percentage: See :class:`BayesianRecommender`.
+        n_restarts: See :class:`BayesianRecommender`.
+        n_raw_samples: See :class:`BayesianRecommender`.
+        max_n_subsets: See :class:`BayesianRecommender`.
+        max_n_subspaces: Deprecated! Use ``max_n_subsets`` instead.
+
+    Returns:
+        BayesianRecommender: Instance of `BayesianRecommender` with provided parameters.
     """
+    from baybe.recommenders.pure.bayesian.base import BayesianRecommender
 
-    max_n_subsets: int = field(default=10, validator=[instance_of(int), ge(1)])
-    """Maximum number of subsets to evaluate when subset-generating constraints are
-    present (e.g., continuous cardinality constraints). If the total number of
-    subsets exceeds this limit, a random subset of that size is sampled for
-    optimization instead of performing an exhaustive search."""
-
-    @property
-    def max_n_subspaces(self) -> int:
-        """Deprecated! Use ``max_n_subsets`` instead."""
+    if max_n_subspaces is not None:
         warnings.warn(
             "'max_n_subspaces' has been renamed to 'max_n_subsets' and will "
             "be removed in a future version.",
             DeprecationWarning,
             stacklevel=2,
         )
-        return self.max_n_subsets
+        max_n_subsets = max_n_subspaces
 
-    @max_n_subspaces.setter
-    def max_n_subspaces(self, value: int) -> None:
-        """Deprecated! Use ``max_n_subsets`` instead."""  # noqa: D401
-        warnings.warn(
-            "'max_n_subspaces' has been renamed to 'max_n_subsets' and will "
-            "be removed in a future version.",
-            DeprecationWarning,
-            stacklevel=2,
-        )
-        self.max_n_subsets = value
-
-    @sampling_percentage.validator
-    def _validate_percentage(  # noqa: DOC101, DOC103
-        self, _: Any, value: float
-    ) -> None:
-        """Validate that the given value is in fact a percentage.
-
-        Raises:
-            ValueError: If ``value`` is not between 0 and 1.
-        """
-        if not 0 <= value <= 1:
-            raise ValueError(
-                f"Hybrid sampling percentage needs to be between 0 and 1 but is {value}"
-            )
+    warnings.warn(
+        "'BotorchRecommender' is deprecated and will be removed in a future version. "
+        "Please use 'BayesianRecommender' instead.",
+        DeprecationWarning,
+        stacklevel=2,
+    )
 
-    @override
-    def __str__(self) -> str:
-        fields = [
-            to_string("Surrogate", self._surrogate_model),
-            to_string(
-                "Acquisition function", self.acquisition_function, single_line=True
+    # TODO: Clean up once more optimizers are implemented.
+    if not sequential_continuous:
+        return BayesianRecommender(
+            surrogate_model=surrogate_model,
+            acquisition_function=acquisition_function,
+            optimizer=GradientOptimizer(
+                n_restarts=n_restarts,
+                n_raw_samples=n_raw_samples,
+                sequential_continuous=sequential_continuous,
             ),
-            to_string("Compatibility", self.compatibility, single_line=True),
-            to_string(
-                "Sequential continuous", self.sequential_continuous, single_line=True
-            ),
-            to_string("Hybrid sampler", self.hybrid_sampler, single_line=True),
-            to_string(
-                "Sampling percentage", self.sampling_percentage, single_line=True
-            ),
-        ]
-        return to_string(self.__class__.__name__, *fields)
-
-    @override
-    def _recommend_discrete(
-        self,
-        subspace_discrete: SubspaceDiscrete,
-        candidates_exp: pd.DataFrame,
-        batch_size: int,
-    ) -> pd.Index:
-        """Generate recommendations from a discrete search space.
-
-        Dispatches to the appropriate optimization routine depending on whether
-        subset constraints are present.
-
-        Args:
-            subspace_discrete: The discrete subspace from which to generate
-                recommendations.
-            candidates_exp: The experimental representation of all discrete candidate
-                points to be considered.
-            batch_size: The size of the recommendation batch.
-
-        Returns:
-            The dataframe indices of the recommended points in the provided
-            experimental representation.
-        """
-        if subspace_discrete.n_subsets > 0:
-            return recommend_discrete_with_subsets(
-                self, subspace_discrete, candidates_exp, batch_size
-            )
-        return recommend_discrete_without_subsets(
-            self, subspace_discrete, candidates_exp, batch_size
+            hybrid_sampler=hybrid_sampler,
+            sampling_percentage=sampling_percentage,
+            n_restarts=n_restarts,
+            n_raw_samples=n_raw_samples,
+            max_n_subsets=max_n_subsets,
         )
-
-    @override
-    def _recommend_continuous(
-        self,
-        subspace_continuous: SubspaceContinuous,
-        batch_size: int,
-    ) -> pd.DataFrame:
-        """Generate recommendations from a continuous search space.
-
-        Args:
-            subspace_continuous: The continuous subspace from which to generate
-                recommendations.
-            batch_size: The size of the recommendation batch.
-
-        Raises:
-            IncompatibleAcquisitionFunctionError: If a non-Monte Carlo acquisition
-                function is used with a batch size > 1.
-
-        Returns:
-            A dataframe containing the recommendations as individual rows.
-        """
-        assert self._objective is not None
-        if (
-            batch_size > 1
-            and not self._get_acquisition_function(self._objective).supports_batching
-        ):
-            raise IncompatibleAcquisitionFunctionError(
-                f"The '{self.__class__.__name__}' only works with Monte Carlo "
-                f"acquisition functions for batch sizes > 1."
-            )
-
-        points, _ = recommend_continuous_torch(self, subspace_continuous, batch_size)
-
-        return pd.DataFrame(points, columns=subspace_continuous.parameter_names)
-
-    @override
-    def _recommend_hybrid(
-        self,
-        searchspace: SearchSpace,
-        candidates_exp: pd.DataFrame,
-        batch_size: int,
-    ) -> pd.DataFrame:
-        """Generate recommendations from a hybrid search space.
-
-        Dispatches to the appropriate optimization routine depending on whether
-        subset constraints are present.
-
-        Args:
-            searchspace: The search space in which the recommendations should be made.
-            candidates_exp: The experimental representation of the candidates
-                of the discrete subspace.
-            batch_size: The size of the calculated batch.
-
-        Returns:
-            The recommended points.
-        """
-        if searchspace.n_subsets > 0:
-            return recommend_hybrid_with_subsets(
-                self, searchspace, candidates_exp, batch_size
-            )
-        return recommend_hybrid_without_subsets(
-            self, searchspace, candidates_exp, batch_size
+    else:
+        return BayesianRecommender(
+            surrogate_model=surrogate_model,
+            acquisition_function=acquisition_function,
+            hybrid_sampler=hybrid_sampler,
+            sampling_percentage=sampling_percentage,
+            n_restarts=n_restarts,
+            n_raw_samples=n_raw_samples,
+            max_n_subsets=max_n_subsets,
         )
 
-    def _optimize_over_subsets(
-        self,
-        subset_callables: Iterable[Callable[[], tuple[Any, Tensor]]],
-    ) -> tuple[Any, Tensor]:
-        """Optimize across subsets and return the result with the best acqf value.
-
-        Each callable performs optimization for one subset configuration and returns
-        a ``(result, acquisition_value)`` tuple. Subsets that raise
-        ``InfeasibilityError`` are silently skipped.
-
-        Args:
-            subset_callables: An iterable of zero-argument callables. Each callable
-                runs the optimization for one subset and returns
-                ``(result, acqf_value)``. It may raise ``InfeasibilityError`` if the
-                subset is infeasible.
-
-        Raises:
-            InfeasibilityError: If none of the subsets has a feasible solution.
-
-        Returns:
-            The result and acquisition value of the best subset.
-        """
-        from botorch.exceptions.errors import InfeasibilityError as BoInfeasibilityError
-
-        results_all: list = []
-        acqf_values_all: list[Tensor] = []
-
-        for optimize_fn in subset_callables:
-            try:
-                result, acqf_value = optimize_fn()
-                results_all.append(result)
-                acqf_values_all.append(acqf_value)
-            except (BoInfeasibilityError, InfeasibilityError):
-                pass
-
-        if not results_all:
-            raise InfeasibilityError(
-                "No feasible solution could be found. Potentially the specified "
-                "constraints are too restrictive, i.e. there may be too many "
-                "constraints or thresholds may have been set too tightly. "
-                "Consider relaxing the constraints to improve the chances "
-                "of finding a feasible solution."
-            )
-
-        best_idx = np.argmax(acqf_values_all)
-        return results_all[best_idx], acqf_values_all[best_idx]
-
 
 # Collect leftover original slotted classes processed by `attrs.define`
 gc.collect()
diff --git a/baybe/recommenders/pure/bayesian/botorch/optimizers/__init__.py b/baybe/recommenders/pure/bayesian/botorch/optimizers/__init__.py
new file mode 100644
index 0000000000..eaf5773900
--- /dev/null
+++ b/baybe/recommenders/pure/bayesian/botorch/optimizers/__init__.py
@@ -0,0 +1,7 @@
+"""Acquisition function optimizers."""
+
+from baybe.recommenders.pure.bayesian.botorch.optimizers.basic import GradientOptimizer
+
+__all__ = [
+    "GradientOptimizer",
+]
diff --git a/baybe/recommenders/pure/bayesian/botorch/optimizers/base.py b/baybe/recommenders/pure/bayesian/botorch/optimizers/base.py
new file mode 100644
index 0000000000..443b716e1f
--- /dev/null
+++ b/baybe/recommenders/pure/bayesian/botorch/optimizers/base.py
@@ -0,0 +1,47 @@
+"""Base protocol for all optimizers."""
+
+from __future__ import annotations
+
+from collections.abc import Callable
+from typing import TYPE_CHECKING, ClassVar, Protocol, TypeAlias, runtime_checkable
+
+from baybe.searchspace import SearchSpace
+from baybe.searchspace.core import SearchSpaceType
+
+if TYPE_CHECKING:
+    from torch import Tensor
+
+    Optimand: TypeAlias = Callable[[Tensor], Tensor]
+    "Type alias for the callable to be optimized."
+
+
+@runtime_checkable
+class OptimizerProtocol(Protocol):
+    """Type protocol specifying the interface optimizers need to implement."""
+
+    # Use slots so that derived classes also remain slotted
+    # See also: https://www.attrs.org/en/stable/glossary.html#term-slotted-classes
+    __slots__ = ()
+
+    compatibility: ClassVar[SearchSpaceType]
+    """Class variable reflecting the search space compatibility."""
+
+    def __call__(
+        self,
+        batch_size: int,
+        acquisition_function: Optimand,
+        searchspace: SearchSpace,
+        fixed_parameters: dict[int, float] | None = None,
+    ) -> tuple[Tensor, Tensor]:
+        """Recommend a batch of points from the given search space.
+
+        Args:
+            batch_size: The size of the recommendation batch.
+            acquisition_function: The acquisition function to be optimized.
+            searchspace: The search space from which to generate recommendations.
+            fixed_parameters: A dictionary mapping parameter indices to fixed values.
+
+        Returns:
+            The recommendations and corresponding acquisition values.
+        """
+        ...
diff --git a/baybe/recommenders/pure/bayesian/botorch/optimizers/basic.py b/baybe/recommenders/pure/bayesian/botorch/optimizers/basic.py
new file mode 100644
index 0000000000..eeae1b69fb
--- /dev/null
+++ b/baybe/recommenders/pure/bayesian/botorch/optimizers/basic.py
@@ -0,0 +1,140 @@
+"""Low-level optimizers of acquisition functions."""
+
+from __future__ import annotations
+
+import gc
+from typing import TYPE_CHECKING, ClassVar, cast
+
+from attrs import define, field, fields
+from attrs.validators import gt, instance_of
+from typing_extensions import override
+
+from baybe.exceptions import IncompatibilityError, IncompatibleSearchSpaceError
+from baybe.recommenders.pure.bayesian.botorch.optimizers.base import OptimizerProtocol
+from baybe.searchspace import SearchSpace
+from baybe.searchspace.core import SearchSpaceType
+from baybe.utils.basic import flatten
+
+if TYPE_CHECKING:
+    from torch import Tensor
+
+    from baybe.recommenders.pure.bayesian.botorch.optimizers.base import Optimand
+
+
+@define(kw_only=True)
+class GradientOptimizer(OptimizerProtocol):
+    """Acquisition function optimizer using gradient-based optimization."""
+
+    # Class variables
+    compatibility: ClassVar[SearchSpaceType] = SearchSpaceType.CONTINUOUS
+    # See base class.
+
+    n_restarts: int = field(validator=[instance_of(int), gt(0)], default=10)
+    """Number of times gradient-based optimization is restarted from different initial
+    points.
+    """
+
+    n_raw_samples: int = field(validator=[instance_of(int), gt(0)], default=64)
+    """Number of raw samples drawn for the initialization heuristic in gradient-based
+    optimization.
+    """
+
+    sequential_continuous: bool = field(default=True)
+    """Flag defining whether to apply sequential greedy or batch optimization in
+    **continuous** search spaces. In discrete/hybrid spaces, sequential greedy
+    optimization is applied automatically.
+    """
+
+    @override
+    def __call__(
+        self,
+        batch_size: int,
+        acquisition_function: Optimand,
+        searchspace: SearchSpace,
+        fixed_parameters: dict[int, float] | None = None,
+    ) -> tuple[Tensor, Tensor]:
+        """Recommend from a search space using gradient-based optimization.
+
+        Args:
+            batch_size: The size of the recommendation batch.
+            acquisition_function: The acquisition function to be optimized.
+            searchspace: The search space from which to generate recommendations.
+            fixed_parameters: A dictionary mapping parameter indices to fixed values.
+
+        Returns:
+            The recommendations and corresponding acquisition values.
+
+        Raises:
+            IncompatibilityError: If the search space has interpoint constraints and the
+                ``sequential_continuous`` flag is set to ``True``.
+            NotImplementedError: If the search space has a discrete component.
+            ValueError: If the search space has cardinality constraints.
+        """
+        import torch
+        from botorch.acquisition import AcquisitionFunction as BoAcquisitionFunction
+        from botorch.optim import optimize_acqf
+
+        if searchspace.type is not self.compatibility:
+            raise IncompatibleSearchSpaceError(
+                f"'{self.__class__.__name__}' currently only supports "
+                f"continuous search spaces."
+            )
+
+        # TODO: Add option for automatic choice once the "settings" PR is merged,
+        #   which ships the necessary machinery
+        if (
+            self.sequential_continuous
+            and searchspace.continuous.has_interpoint_constraints
+        ):
+            raise IncompatibilityError(
+                f"Setting the "
+                f"'{fields(self.__class__).sequential_continuous.name}' "
+                f"flag to ``True`` while interpoint constraints are present in the "
+                f"continuous subspace is not supported. "
+            )
+
+        if not searchspace.discrete.is_empty:
+            raise NotImplementedError(
+                "Gradient-based optimization is not implemented "
+                "for non-empty discrete search spaces."
+            )
+
+        if searchspace.continuous.n_subsets > 0:
+            raise ValueError(
+                f"'{self.__class__.__name__}' "
+                f"expects a continuous subspace without cardinality constraints."
+            )
+
+        points, acqf_values = optimize_acqf(
+            acq_function=cast(BoAcquisitionFunction, acquisition_function),
+            bounds=torch.from_numpy(
+                searchspace.continuous.comp_rep_bounds.to_numpy(copy=True)
+            ),
+            q=batch_size,
+            num_restarts=self.n_restarts,
+            raw_samples=self.n_raw_samples,
+            fixed_features=fixed_parameters or None,
+            equality_constraints=flatten(
+                c.to_botorch(
+                    searchspace.continuous.parameters,
+                    batch_size=batch_size if c.is_interpoint else None,
+                )
+                for c in searchspace.continuous.constraints_lin_eq
+            )
+            or None,
+            inequality_constraints=flatten(
+                c.to_botorch(
+                    searchspace.continuous.parameters,
+                    batch_size=batch_size if c.is_interpoint else None,
+                )
+                for c in searchspace.continuous.constraints_lin_ineq
+            )
+            or None,
+            sequential=self.sequential_continuous,
+        )
+
+        return points, acqf_values
+
+
+# Collect leftover original slotted classes processed by `attrs.define`
+gc.collect()
diff --git a/baybe/recommenders/pure/bayesian/botorch/continuous.py b/baybe/recommenders/pure/bayesian/continuous.py
similarity index 73%
rename from baybe/recommenders/pure/bayesian/botorch/continuous.py
rename to baybe/recommenders/pure/bayesian/continuous.py
index 4e43a247fb..42cb8d9b90 100644
--- a/baybe/recommenders/pure/bayesian/botorch/continuous.py
+++ b/baybe/recommenders/pure/bayesian/continuous.py
@@ -1,4 +1,4 @@
-"""Continuous recommendation routines for BotorchRecommender."""
+"""Continuous recommendation routines for BayesianRecommender."""
 
 from __future__ import annotations
 
@@ -7,25 +7,21 @@
 from typing import TYPE_CHECKING
 
 import pandas as pd
-from attrs import fields
 
 from baybe.constraints.utils import is_cardinality_fulfilled
-from baybe.exceptions import (
-    IncompatibilityError,
-    MinimumCardinalityViolatedWarning,
-)
+from baybe.exceptions import MinimumCardinalityViolatedWarning
 from baybe.parameters.numerical import _FixedNumericalContinuousParameter
 from baybe.searchspace import SubspaceContinuous
-from baybe.utils.basic import flatten
+from baybe.searchspace.core import SearchSpace
 
 if TYPE_CHECKING:
     from torch import Tensor
 
-    from baybe.recommenders.pure.bayesian.botorch.core import BotorchRecommender
+    from baybe.recommenders.pure.bayesian.base import BayesianRecommender
 
 
 def recommend_continuous_torch(
-    recommender: BotorchRecommender,
+    recommender: BayesianRecommender,
     subspace_continuous: SubspaceContinuous,
     batch_size: int,
 ) -> tuple[Tensor, Tensor]:
@@ -41,7 +37,7 @@ def recommend_continuous_torch(
 
 
 def recommend_continuous_with_cardinality_constraints(
-    recommender: BotorchRecommender,
+    recommender: BayesianRecommender,
     subspace_continuous: SubspaceContinuous,
     batch_size: int,
 ) -> tuple[Tensor, Tensor]:
@@ -129,7 +125,7 @@ def optimize() -> tuple[Tensor, Tensor]:
 
 
 def recommend_continuous_without_cardinality_constraints(
-    recommender: BotorchRecommender,
+    recommender: BayesianRecommender,
     subspace_continuous: SubspaceContinuous,
     batch_size: int,
 ) -> tuple[Tensor, Tensor]:
@@ -147,9 +143,6 @@ def recommend_continuous_without_cardinality_constraints(
     Raises:
         ValueError: If the continuous search space has cardinality constraints.
     """
-    import torch
-    from botorch.optim import optimize_acqf
-
     if subspace_continuous.n_subsets > 0:
         raise ValueError(
             f"'{recommend_continuous_without_cardinality_constraints.__name__}' "
@@ -162,51 +155,14 @@ def recommend_continuous_without_cardinality_constraints(
         if isinstance(p, _FixedNumericalContinuousParameter)
     }
 
-    # TODO: Add option for automatic choice once the "settings" PR is merged,
-    #   which ships the necessary machinery
-    if (
-        recommender.sequential_continuous
-        and subspace_continuous.has_interpoint_constraints
-    ):
-        from baybe.recommenders.pure.bayesian.botorch.core import BotorchRecommender
-
-        raise IncompatibilityError(
-            f"Setting the "
-            f"'{fields(BotorchRecommender).sequential_continuous.name}' "
-            f"flag to ``True`` while interpoint constraints are present in the "
-            f"continuous subspace is not supported. "
-        )
-
     # NOTE: The explicit `or None` conversion is added as an additional safety net
     #   because it is unclear if the corresponding presence checks for these
     #   arguments is correctly implemented in all invoked BoTorch subroutines.
     #   For details: https://github.com/pytorch/botorch/issues/2042
-    points, acqf_values = optimize_acqf(
-        acq_function=recommender._botorch_acqf,
-        bounds=torch.from_numpy(
-            subspace_continuous.comp_rep_bounds.to_numpy(copy=True)
-        ),
-        q=batch_size,
-        num_restarts=recommender.n_restarts,
-        raw_samples=recommender.n_raw_samples,
-        fixed_features=fixed_parameters or None,
-        equality_constraints=flatten(
-            c.to_botorch(
-                subspace_continuous.parameters,
-                batch_size=batch_size if c.is_interpoint else None,
-            )
-            for c in subspace_continuous.constraints_lin_eq
-        )
-        or None,
-        inequality_constraints=flatten(
-            c.to_botorch(
-                subspace_continuous.parameters,
-                batch_size=batch_size if c.is_interpoint else None,
-            )
-            for c in subspace_continuous.constraints_lin_ineq
-        )
-        or None,
-        sequential=recommender.sequential_continuous,
+    points, acqf_values = recommender.optimizer(
+        batch_size=batch_size,
+        acquisition_function=recommender._botorch_acqf,
+        searchspace=SearchSpace(continuous=subspace_continuous),
+        fixed_parameters=fixed_parameters,
     )
-    assert acqf_values is not None  # for mypy; guaranteed by optimize_acqf defaults
     return points, acqf_values
diff --git a/baybe/recommenders/pure/bayesian/botorch/discrete.py b/baybe/recommenders/pure/bayesian/discrete.py
similarity index 95%
rename from baybe/recommenders/pure/bayesian/botorch/discrete.py
rename to baybe/recommenders/pure/bayesian/discrete.py
index a5f92d04d0..9d970d5585 100644
--- a/baybe/recommenders/pure/bayesian/botorch/discrete.py
+++ b/baybe/recommenders/pure/bayesian/discrete.py
@@ -1,4 +1,4 @@
-"""Discrete recommendation routines for BotorchRecommender."""
+"""Discrete recommendation routines for BayesianRecommender."""
 
 from __future__ import annotations
 
@@ -15,11 +15,11 @@
 if TYPE_CHECKING:
     from torch import Tensor
 
-    from baybe.recommenders.pure.bayesian.botorch.core import BotorchRecommender
+    from baybe.recommenders.pure.bayesian.base import BayesianRecommender
 
 
 def recommend_discrete_with_subsets(
-    recommender: BotorchRecommender,
+    recommender: BayesianRecommender,
     subspace_discrete: SubspaceDiscrete,
     candidates_exp: pd.DataFrame,
     batch_size: int,
@@ -76,7 +76,7 @@ def optimize() -> tuple[pd.Index, Tensor]:
 
 
 def recommend_discrete_without_subsets(
-    recommender: BotorchRecommender,
+    recommender: BayesianRecommender,
     subspace_discrete: SubspaceDiscrete,
     candidates_exp: pd.DataFrame,
     batch_size: int,
diff --git a/baybe/recommenders/pure/bayesian/botorch/hybrid.py b/baybe/recommenders/pure/bayesian/hybrid.py
similarity index 97%
rename from baybe/recommenders/pure/bayesian/botorch/hybrid.py
rename to baybe/recommenders/pure/bayesian/hybrid.py
index 0b017a92dd..c82babeea9 100644
--- a/baybe/recommenders/pure/bayesian/botorch/hybrid.py
+++ b/baybe/recommenders/pure/bayesian/hybrid.py
@@ -1,4 +1,4 @@
-"""Hybrid recommendation routines for BotorchRecommender."""
+"""Hybrid recommendation routines for BayesianRecommender."""
 
 from __future__ import annotations
 
@@ -24,11 +24,11 @@
 if TYPE_CHECKING:
     from torch import Tensor
 
-    from baybe.recommenders.pure.bayesian.botorch.core import BotorchRecommender
+    from baybe.recommenders.pure.bayesian.base import BayesianRecommender
 
 
 def recommend_hybrid_without_subsets(
-    recommender: BotorchRecommender,
+    recommender: BayesianRecommender,
     searchspace: SearchSpace,
     candidates_exp: pd.DataFrame,
     batch_size: int,
@@ -162,7 +162,7 @@ def recommend_hybrid_without_subsets(
 
 
 def recommend_hybrid_with_subsets(
-    recommender: BotorchRecommender,
+    recommender: BayesianRecommender,
     searchspace: SearchSpace,
     candidates_exp: pd.DataFrame,
     batch_size: int,
diff --git a/examples/Basics/recommenders.py b/examples/Basics/recommenders.py
index 14a87731d1..39309022e3 100644
--- a/examples/Basics/recommenders.py
+++ b/examples/Basics/recommenders.py
@@ -20,7 +20,7 @@
 from baybe.objectives import SingleTargetObjective
 from baybe.parameters import NumericalDiscreteParameter, SubstanceParameter
 from baybe.recommenders import (
-    BotorchRecommender,
+    BayesianRecommender,
     RandomRecommender,
     TwoPhaseMetaRecommender,
 )
@@ -96,7 +96,7 @@
 
 recommender = TwoPhaseMetaRecommender(
     initial_recommender=INITIAL_RECOMMENDER,
-    recommender=BotorchRecommender(
+    recommender=BayesianRecommender(
         surrogate_model=SURROGATE_MODEL, acquisition_function=ACQ_FUNCTION
     ),
 )
diff --git a/examples/Constraints_Continuous/interpoint.py b/examples/Constraints_Continuous/interpoint.py
index 0bfe1584f6..30f01a9d79 100644
--- a/examples/Constraints_Continuous/interpoint.py
+++ b/examples/Constraints_Continuous/interpoint.py
@@ -30,7 +30,7 @@
 from baybe import Campaign, active_settings
 from baybe.constraints import ContinuousLinearConstraint
 from baybe.parameters import NumericalContinuousParameter
-from baybe.recommenders import BotorchRecommender
+from baybe.recommenders import BayesianRecommender
 from baybe.recommenders.meta.sequential import TwoPhaseMetaRecommender
 from baybe.searchspace import SearchSpace
 from baybe.targets import NumericalTarget
@@ -122,7 +122,7 @@
 )
 objective = NumericalTarget(name="Reaction_Yield").to_objective()
 recommender = TwoPhaseMetaRecommender(
-    recommender=BotorchRecommender(sequential_continuous=False)
+    recommender=BayesianRecommender(sequential_continuous=False)
 )
 campaign = Campaign(
     searchspace=searchspace,
diff --git a/examples/Custom_Hooks/campaign_stopping.py b/examples/Custom_Hooks/campaign_stopping.py
index 208cbe60f9..28339f5317 100644
--- a/examples/Custom_Hooks/campaign_stopping.py
+++ b/examples/Custom_Hooks/campaign_stopping.py
@@ -26,7 +26,7 @@
 from baybe.objectives.base import Objective
 from baybe.parameters import NumericalDiscreteParameter, SubstanceParameter
 from baybe.recommenders import (
-    BotorchRecommender,
+    BayesianRecommender,
     RandomRecommender,
     TwoPhaseMetaRecommender,
 )
@@ -87,7 +87,7 @@
 objective = SingleTargetObjective(target=NumericalTarget(name="yield"))
 
 recommender = TwoPhaseMetaRecommender(
-    initial_recommender=RandomRecommender(), recommender=BotorchRecommender()
+    initial_recommender=RandomRecommender(), recommender=BayesianRecommender()
 )
 
 ### Simulating the Uninterrupted Campaigns
@@ -126,7 +126,7 @@ class CampaignStoppedException(Exception):
 
 
 def stop_on_PI(
-    self: BotorchRecommender,
+    self: BayesianRecommender,
     searchspace: SearchSpace,
     objective: Objective | None = None,
     measurements: pd.DataFrame | None = None,
@@ -155,8 +155,8 @@ def stop_on_PI(
 
 # Now, we attach the hook to the ``recommend`` function of our recommender class:
 
-BotorchRecommender.recommend = register_hooks(
-    BotorchRecommender.recommend, post_hooks=[stop_on_PI]
+BayesianRecommender.recommend = register_hooks(
+    BayesianRecommender.recommend, post_hooks=[stop_on_PI]
 )
 active_settings.parallelize_simulation_runs = False
 
@@ -178,7 +178,7 @@ def stop_on_PI(
 # hook and assign it to a fresh copy of the campaign:
 
 recommender_with_hook = TwoPhaseMetaRecommender(
-    initial_recommender=RandomRecommender(), recommender=BotorchRecommender()
+    initial_recommender=RandomRecommender(), recommender=BayesianRecommender()
 )
 campaign_with_hook = Campaign(searchspace, objective, recommender)
 
diff --git a/examples/Custom_Hooks/probability_of_improvement.py b/examples/Custom_Hooks/probability_of_improvement.py
index 70511fffeb..49b7534c37 100644
--- a/examples/Custom_Hooks/probability_of_improvement.py
+++ b/examples/Custom_Hooks/probability_of_improvement.py
@@ -4,7 +4,7 @@
 # {func}`register_hooks <baybe.utils.basic.register_hooks>` utility can be used to
 # extract the *Probability of Improvement (PI)* from a running campaign:
 # * We define a hook that is compatible with the
-#   {meth}`BotorchRecommender.recommend <baybe.recommenders.pure.bayesian.botorch.core.BotorchRecommender.recommend>`
+#   {meth}`BayesianRecommender.recommend <baybe.recommenders.pure.bayesian.botorch.core.BayesianRecommender.recommend>`
 #   interface and lets us extract the PI achieved after each experimental iteration,
 # * attach the hook to the recommender driving our campaign,
 # * and plot the evolving PI values after campaign completion.
@@ -29,7 +29,7 @@
 from baybe.objectives.base import Objective
 from baybe.parameters import NumericalDiscreteParameter
 from baybe.recommenders import (
-    BotorchRecommender,
+    BayesianRecommender,
     RandomRecommender,
     TwoPhaseMetaRecommender,
 )
@@ -67,7 +67,7 @@
 
 
 def extract_pi(
-    self: BotorchRecommender,
+    self: BayesianRecommender,
     searchspace: SearchSpace,
     objective: Objective | None = None,
     measurements: pd.DataFrame | None = None,
@@ -91,12 +91,12 @@ def extract_pi(
 
 # Next, we create our recommender and monkeypatch its `recommend` method:
 
-bayesian_recommender = BotorchRecommender(
+bayesian_recommender = BayesianRecommender(
     surrogate_model=GaussianProcessSurrogate(),
 )
 bayesian_recommender.recommend = MethodType(
     register_hooks(
-        BotorchRecommender.recommend,
+        BayesianRecommender.recommend,
         post_hooks=[extract_pi],
     ),
     bayesian_recommender,
@@ -107,7 +107,7 @@ def extract_pi(
 )
 
 # In this example, we use `MethodType` to bind the
-# {meth}`BotorchRecommender.recommend <baybe.recommenders.pure.bayesian.botorch.core.BotorchRecommender.recommend>`
+# {meth}`BayesianRecommender.recommend <baybe.recommenders.pure.bayesian.botorch.core.BayesianRecommender.recommend>`
 # **function** with our hook.
 # For more information, we refer to the [`basic example`](./basics.md) explaining the
 # hook mechanics.
diff --git a/examples/Custom_Surrogates/custom_pretrained.py b/examples/Custom_Surrogates/custom_pretrained.py
index 5e408fd4b4..b6dc3bf46a 100644
--- a/examples/Custom_Surrogates/custom_pretrained.py
+++ b/examples/Custom_Surrogates/custom_pretrained.py
@@ -19,7 +19,7 @@
 from baybe.objectives import SingleTargetObjective
 from baybe.parameters import NumericalDiscreteParameter
 from baybe.recommenders import (
-    BotorchRecommender,
+    BayesianRecommender,
     FPSRecommender,
     TwoPhaseMetaRecommender,
 )
@@ -102,7 +102,7 @@
     searchspace=SearchSpace.from_product(parameters=parameters, constraints=None),
     objective=SingleTargetObjective(target=NumericalTarget(name="Yield")),
     recommender=TwoPhaseMetaRecommender(
-        recommender=BotorchRecommender(surrogate_model=surrogate_model),
+        recommender=BayesianRecommender(surrogate_model=surrogate_model),
         initial_recommender=FPSRecommender(),
     ),
 )
diff --git a/examples/Custom_Surrogates/surrogate_params.py b/examples/Custom_Surrogates/surrogate_params.py
index d9d4852fbe..b890f24d4c 100644
--- a/examples/Custom_Surrogates/surrogate_params.py
+++ b/examples/Custom_Surrogates/surrogate_params.py
@@ -19,7 +19,7 @@
     SubstanceParameter,
 )
 from baybe.recommenders import (
-    BotorchRecommender,
+    BayesianRecommender,
     FPSRecommender,
     TwoPhaseMetaRecommender,
 )
@@ -84,7 +84,7 @@
     searchspace=SearchSpace.from_product(parameters=parameters, constraints=None),
     objective=SingleTargetObjective(target=NumericalTarget(name="Yield")),
     recommender=TwoPhaseMetaRecommender(
-        recommender=BotorchRecommender(surrogate_model=surrogate_model),
+        recommender=BayesianRecommender(surrogate_model=surrogate_model),
         initial_recommender=FPSRecommender(),
     ),
 )
diff --git a/examples/Multi_Armed_Bandit/bernoulli_multi_armed_bandit.py b/examples/Multi_Armed_Bandit/bernoulli_multi_armed_bandit.py
index d8c80317cb..9728952b09 100644
--- a/examples/Multi_Armed_Bandit/bernoulli_multi_armed_bandit.py
+++ b/examples/Multi_Armed_Bandit/bernoulli_multi_armed_bandit.py
@@ -19,7 +19,7 @@
 from baybe.acquisition.base import AcquisitionFunction
 from baybe.parameters import CategoricalParameter
 from baybe.recommenders import (
-    BotorchRecommender,
+    BayesianRecommender,
     RandomRecommender,
     TwoPhaseMetaRecommender,
 )
@@ -124,7 +124,7 @@ def simulate(acqf: AcquisitionFunction) -> SimulationResult:
     """Simulate the campaign with with the given acquisition function."""
     recommender = TwoPhaseMetaRecommender(
         initial_recommender=RandomRecommender(),
-        recommender=BotorchRecommender(
+        recommender=BayesianRecommender(
             surrogate_model=surrogate, acquisition_function=acqf
         ),
     )
diff --git a/examples/Searchspaces/hybrid_space.py b/examples/Searchspaces/hybrid_space.py
index 7d2b82d832..3013891d91 100644
--- a/examples/Searchspaces/hybrid_space.py
+++ b/examples/Searchspaces/hybrid_space.py
@@ -114,7 +114,7 @@
 # Here, we explicitly create a recommender object to use the `NaiveHybridSpaceRecommender`.
 # The keywords `disc_recommender` and `cont_recommender` can be used to select different
 # recommenders for the corresponding subspaces.
-# We use the default choices, which is the `BotorchRecommender`.
+# We use the default choices, which is the `BayesianRecommender`.
 
 hybrid_recommender = TwoPhaseMetaRecommender(recommender=NaiveHybridSpaceRecommender())
 
diff --git a/examples/Serialization/basic_serialization.py b/examples/Serialization/basic_serialization.py
index 77f48dac79..001853f202 100644
--- a/examples/Serialization/basic_serialization.py
+++ b/examples/Serialization/basic_serialization.py
@@ -17,7 +17,7 @@
     NumericalDiscreteParameter,
 )
 from baybe.recommenders import (
-    BotorchRecommender,
+    BayesianRecommender,
     FPSRecommender,
     TwoPhaseMetaRecommender,
 )
@@ -49,7 +49,7 @@
     searchspace=SearchSpace.from_product(parameters=parameters, constraints=None),
     objective=SingleTargetObjective(target=NumericalTarget(name="Yield")),
     recommender=TwoPhaseMetaRecommender(
-        recommender=BotorchRecommender(),
+        recommender=BayesianRecommender(),
         initial_recommender=FPSRecommender(),
     ),
 )
diff --git a/examples/Serialization/create_from_config.py b/examples/Serialization/create_from_config.py
index 4385599090..9ab67df71c 100644
--- a/examples/Serialization/create_from_config.py
+++ b/examples/Serialization/create_from_config.py
@@ -72,7 +72,7 @@
             "type": "FPSRecommender"
         },
         "recommender": {
-            "type": "BotorchRecommender",
+            "type": "BayesianRecommender",
             "surrogate_model": {
                 "type": "GaussianProcessSurrogate"
             },
diff --git a/examples/Serialization/validate_config.py b/examples/Serialization/validate_config.py
index 166266acf1..3377e28213 100644
--- a/examples/Serialization/validate_config.py
+++ b/examples/Serialization/validate_config.py
@@ -71,7 +71,7 @@
             "type": "FPSRecommender"
         },
         "recommender": {
-            "type": "BotorchRecommender",
+            "type": "BayesianRecommender",
             "surrogate_model": {
                 "type": "GaussianProcessSurrogate"
             },
@@ -140,7 +140,7 @@
             "type": "FPSRecommender"
         },
         "recommender": {
-            "type": "BotorchRecommender",
+            "type": "BayesianRecommender",
             "surrogate_model": {
                 "type": "GaussianProcessSurrogate"
             },
diff --git a/streamlit/surrogate_models.py b/streamlit/surrogate_models.py
index 3630f1d2f8..f7e3ebcf31 100644
--- a/streamlit/surrogate_models.py
+++ b/streamlit/surrogate_models.py
@@ -20,7 +20,7 @@
 from baybe.acquisition.base import AcquisitionFunction
 from baybe.exceptions import IncompatibleSurrogateError
 from baybe.parameters import NumericalDiscreteParameter, TaskParameter
-from baybe.recommenders import BotorchRecommender
+from baybe.recommenders import BayesianRecommender
 from baybe.searchspace import SearchSpace
 from baybe.surrogates import CustomONNXSurrogate, GaussianProcessSurrogate
 from baybe.surrogates.base import Surrogate
@@ -272,7 +272,7 @@ def make_surrogate():
             task_meas = measurements[measurements["task"] == task_name][
                 ["x", "y"]
             ].reset_index(drop=True)
-            task_recommender = BotorchRecommender(
+            task_recommender = BayesianRecommender(
                 surrogate_model=make_surrogate(),
                 acquisition_function=acqf,
             )
@@ -296,7 +296,7 @@ def make_surrogate():
             )
     else:
         # Single recommender (single task, or multi-task with transfer learning)
-        recommender = BotorchRecommender(
+        recommender = BayesianRecommender(
             surrogate_model=make_surrogate(),
             acquisition_function=acqf,
         )