What toolbox enables flexible Bayesian optimization in R?

A modular toolbox for building and running BO

Bayesian optimization (BO) targets objective functions whose internal structure is unknown and whose gradients are unavailable, so the only way to learn about them is to evaluate inputs and observe outputs. The paper presents mlr3mbo as an extensible R toolbox for this setting, built from object-oriented components that can be combined and extended. It sits inside the mlr3 ecosystem and builds on bbotk for black-box objectives and optimization and on paradox for defining parameter spaces.

The toolbox serves three distinct user groups. Researchers and students can assemble custom BO variants from the building blocks; practitioners can run a standard but highly configurable BO variant on black-box problems; and hyperparameter optimization experiments can use mlr3mbo as a BO backend through its integration into mlr3 and mlr3tuning. The advanced features listed by the authors include arbitrary regression models as surrogate models, handling of mixed and hierarchical search spaces, parallel multi-point proposals in batch and asynchronous modes, multi-objective optimization, and robust error handling.

Surrogate and acquisition choices follow search-space structure

BO approximates the unknown objective with a probabilistic surrogate model and then uses an acquisition function to choose the next point to evaluate. The surrogate is a cheap stand-in for the expensive objective, fitted to the data observed so far. The paper states that for purely numeric search spaces, Gaussian processes are usually the recommended choice because they provide posterior uncertainty and can be adapted through a covariance kernel. When the space includes categorical parameters or hierarchical dependencies, random forests are a viable alternative because they handle such parameters directly, though their uncertainty must be estimated through approaches such as the ensemble standard deviation or the law of total variance.

The acquisition function converts the surrogate's posterior mean and standard deviation into a single criterion that balances exploitation of low predicted loss against exploration of uncertain regions. Expected improvement, the most commonly used acquisition function, estimates the expected gain over the current best observation with a lower bound at zero. The paper also describes the lower confidence bound, which subtracts a scaled posterior standard deviation from the posterior mean, and LogEI, which modifies expected improvement for surrogates trained on log-transformed targets. New points are then obtained by optimizing the acquisition function, which is usually highly multi-modal; the paper lists random search, DIRECT, restarted local search, CMA-ES, and L-BFGS-B as common approaches.

Coordinate descent produced separate defaults for two regimes

Because BO has different requirements for purely numeric and mixed-hierarchical problems, the authors defined separate configuration search spaces for the two regimes. They discretized the numeric parameters of mlr3mbo into a small set of representative values so that each search space became purely categorical, then ran a greedy coordinate-wise local search that changes one parameter at a time and accepts the best configuration in each one-parameter exchange neighborhood. Starting from a simple or standard configuration meant that more complex settings were added only when empirically justified.

Each configuration was evaluated on the YAHPO-SO v1 suite from YAHPO Gym, a surrogate-based benchmark of hyperparameter optimization tasks with varying metrics, dimensionalities, and complexities. Each instance was repeated 30 times with different random seeds, using a budget of ⌈100 + 40√d⌉ evaluations for d dimensions. Performance was scored with the random search normalized score (RSNS), where 0 corresponds to random search with that same budget and 1 corresponds to random search with 1,000,000 evaluations; the score was computed from mean performance across the 30 runs and then averaged over instances.

Tuned defaults perform on par with established optimizers

The paper benchmarks mlr3mbo against a wide range of established optimizers, including HEBO, SMAC3, Ax, and Optuna, and reports that it performs on par with the state of the art on the numeric and mixed-hierarchical YAHPO Gym benchmarks while offering good runtime behavior. This is a comparison of aggregate benchmark performance, not a claim that mlr3mbo wins on every instance or every dimension.

The discussion identifies where the default configurations diverge from competitors. On numeric instances, mlr3mbo's default uses a Gaussian process surrogate like Ax and SMAC4BB, but relies on the lower confidence bound acquisition function where SMAC3 and Ax use expected improvement. On mixed-hierarchical instances, both SMAC4HPO and mlr3mbo use a random forest surrogate and local search for acquisition optimization, but SMAC3 uses only 10 trees and estimates variance with the ensemble standard deviation, whereas mlr3mbo uses 500 trees and the law of total variance. When mlr3mbo was configured with only 10 trees, it performed notably worse, which the authors read as evidence that very small forests give insufficiently stable uncertainty estimates for effective acquisition optimization.

Acquisition optimizer strength matters more than exact optimization

The authors report that strong acquisition function optimizers with sufficiently large budgets matter, and that CMA-ES and local search are the best choices for numeric and mixed-hierarchical instances respectively. They also note that log-transformed targets help. These findings are presented as consistent with prior studies cited in the paper.

The paper records a counterintuitive result about gradient-based optimization. mlr3mbo performed much worse when using L-BFGS-B as the acquisition optimizer, even though Ax successfully relies on L-BFGS-B. The authors suggest this may be because Ax uses automatic differentiation to compute acquisition gradients, whereas mlr3mbo estimates them via finite differences. They add that L-BFGS-B was a strong candidate in the first coordinate descent iteration, indicating that parameters chosen in later iterations interact poorly with L-BFGS-B and that the final configuration is better tuned to CMA-ES and local search. The paper also notes that mlr3mbo benefits from a simple random initial design, whereas Ax and SMAC4BB rely on space-filling designs.

What the benchmark comparison does and does not establish

The coordinate descent search and the performance comparison were conducted on the same set of benchmark instances. The authors acknowledge that this could in principle overfit the mlr3mbo configuration to the YAHPO instances, but argue that the suite is large and diverse, covering AutoML pipelines and a variety of ML algorithms across different dimensionalities, which they believe substantially reduces that risk. They also note that this practice is consistent with how HEBO, SMAC3, BoTorch/Ax, and Optuna are developed and benchmarked on the same families of problems rather than strictly held-out external families.

The coordinate descent procedure itself has two stated limitations. Because it optimizes one parameter at a time while holding the others fixed, it can produce suboptimal configurations when parameter interactions are strong; for example, improving the surrogate to obtain a better uncertainty estimate only helps if the acquisition function exploits that uncertainty. The final configuration is also influenced by the starting configuration. In addition, the RSNS metric is noisy, so the procedure may move in random directions in the last iterations despite averaging over many instances and 30 repetitions; the authors handled this by running the search long enough and fixing a stopping point by inspection rather than through a formal sequential test.

Match the surrogate and acquisition optimizer to the search space

The paper's most decision-relevant split is between numeric and mixed-hierarchical search spaces. For numeric instances, the evidence points to a Gaussian process surrogate with the lower confidence bound acquisition function and CMA-ES as the acquisition optimizer; for mixed-hierarchical instances, it points to a random forest with the law of total variance uncertainty estimator and local search. The 10-tree versus 500-tree comparison makes forest size a concrete setting to check rather than an implementation detail, because the smaller forest performed notably worse in the authors' configuration.

For readers deciding whether to adopt mlr3mbo, the relevant claim is parity rather than superiority: the tuned defaults matched established optimizers on the YAHPO Gym benchmarks, and the toolbox's distinguishing value lies in the ability to construct custom BO variants from its building blocks. The defaults are a starting point validated on that benchmark suite, and the authors' own account of coordinate descent means the configuration is shaped by the chosen start point and by one-parameter-at-a-time search.

Sources used in this answer