arXiv CSopen access

Frugal Bayesian Optimization: Scalable Surrogates for Data- and Resource-Limited Discovery

Frugal Bayesian Optimization: Scalable Surrogates for Data- and Resource-Limited Discovery Panagiotis Krokidas1*, Christoforos Rekatsinas1,2 , Vassilis Sioros1,3 , Grigorios M. Chatziathanasiou4 , Efi-Maria Papia5,6 , George Giannakopoulos1,7

arXiv:2607.29225v1 [cs.LG] 31 Jul 2026

1*

Institute of Informatics and Telecommunications, National Centre for Scientific Research "Demokritos", Agia Paraskevi, Greece. 2 Department of Mechanical Engineering and Aeronautics, University of Patras, Patras Greece. 3 Department of Informatics and Telecommunications, National and Kapodistrian University of Athens, Athens, Greece. 4 School of Mechanical Engineering, Hellenic Mediterranean University, Heraklion, Crete, Greece. 5* Institute of Nanoscience and Nanotechnology, National Centre for Scientific Research "Demokritos", Agia Paraskevi, Greece. 6 Department of Physics, National and Kapodistrian University of Athens, Athens, Greece. 7 SciFY PNPC, Agia Paraskevi, Greece.

*Corresponding author(s). E-mail(s): [email protected]; Contributing authors: [email protected]; [email protected]; [email protected]; [email protected]; [email protected]; Abstract Bayesian Optimization (BO) is widely adopted for data-efficient optimization in scientific and engineering applications, yet its computational cost is rarely evaluated alongside optimization performance. Here we present a systematic, compute-aware study of BO that redefines surrogate evaluation along two axes: optimization quality and computational frugality. Across eight benchmark functions and nine real-world datasets spanning materials science, mechanics, robotics, chemistry, and machine learning, we benchmark four surrogate models—Gaussian Processes, Random Forests, NGBoost, and Bayesian Adaptive Spline Surfaces.

1

We show that Gaussian Process–based BO consistently incurs the highest time and memory overhead without delivering superior optimization or sample efficiency. In contrast, scalable alternatives achieve equal or better performance at a fraction of the computational cost. Motivated by these findings, we introduce a surrogate-recommendation framework that predicts the most suitable BO surrogate from inexpensive dataset characteristics. Together, these results establish FruBO as a reproducible, compute-aware baseline for Bayesian Optimization and provide practical guidance for surrogate selection under limited computational and experimental budgets. Keywords: Bayesian optimization, surrogate modeling, data efficiency, computational scalability, resource-limited discovery

With the advent of machine learning (ML), data-driven methods have become integral to daily life. Production systems (recommendation engines, spam filters, driverassistance) thrive on abundant data, enabling models to interpolate within well-covered distributions. Scientific discovery faces a different regime: data are scarce, and the goal is not to model the average case but to extrapolate toward rare, exceptional solutions. Yet extrapolation is notoriously brittle for models trained on random samples that under-cover the regions of highest interest. To move beyond passive sampling, the community has revisited active learning (AL) [1]: rather than consuming data indiscriminately, the model proposes which sample to acquire next (e.g., a new experiment), guided by informativeness criteria. While AL can yield good predictors with fewer measurements, its iterative loop is typically geared toward representing the space broadly. This does not guarantee accuracy precisely where scientists care most: near the global optimum or within a top-N subset of highvalue candidates. In practice, the scientist’s question is rarely “what happens across the entire design space?” but rather “where is the best solution in the unexplored design space?” This shifts the problem from broad representation to targeted optimization. Bayesian Optimization (BO) [2] addresses this optimization-centric goal by coupling a surrogate model with an acquisition rule to navigate design spaces efficiently. In materials and chemistry, BO has steered searches toward high-performing regions, e.g., for material and drug discovery [3–6]. However, standard BO pipelines rely on Gaussian processes (GPs), whose per-iteration updates involve covariance factorizations with O(N 3 ) time and memory that grows with N (often effectively O(N 2 ) due to kernel matrices) [7, 8]. As acquisitions accumulate, wall-clock time and memory can become prohibitive on modest hardware, pushing users toward specialized GPUs or HPC clusters. Although BO often outperforms state-of-the-art non-BO optimizers, such as evolutionary strategies (e.g., genetic algorithm, CMA-ES) and particle swarm optimization, in terms of sample efficiency at small evaluation budgets [9–11], the GP overhead can dominate as N grows. This is especially true when objective evaluations are inexpensive or highly parallelizable, turning surrogate updates into the bottleneck and limiting BO’s practical use over other optimizers [10, 11]. A large body of work has sought to mitigate the scaling bottlenecks of GP-based BO—via bounded training sets and zooming strategies that progressively narrow the search space [12], focalized sparse-GP objectives that concentrate learning on high-acquisition regions [13], goal-aware acquisition functions that optimize for userdefined subsets rather than the global optimum [3], fixed-memory buffer architectures 2

that cap kernel growth [14], and even hardware acceleration through in-memory or neuromorphic computing [15]. These approaches can substantially delay the onset of cubic cost, but once we commit to GP posteriors we inevitably inherit additional machinery—sparse approximations, variational bounds, or posterior sampling—that increases both methodological and engineering complexity. Such layers of approximation often require expert tuning and careful calibration of hyperparameters, while still failing to eliminate the fundamental memory growth of GP kernels. This motivates exploring surrogates that render BO’s cost nearly O(N ) in time with respect to the number of acquisitions N , while keeping memory usage effectively constant (bounded) in N . Recent studies have indeed begun to explore non-GP surrogates. Gradient-informed Bayesian neural networks trained with stochastic-gradient MCMC can incorporate derivative information and improve sample efficiency, yet explicit wall-clock and memory comparisons against GPs are rarely quantified [16]. The only compute evidence provided is a report of essentially identical training times across three functions—McCormick (2D), Rosenbrock (4D), and Hartmann–6—of about 23 s per training run, indicating stability with respect to dimensionality but no scaling study versus dataset size, nor memory measurements or head-to-head comparisons with GPs. Broader benchmarks examining finite- and infinite-width BNNs, Laplace approximations, deep kernel learning (DKL), and ensemble models show that GPs are not uniquely superior for identifying optima; performance is problem-dependent, with HMC-based BNNs excelling in some regimes and DKL remaining consistently competitive [17]. However, these studies prioritize optimization outcomes and do not profile wall-clock time or memory usage. Lakshminarayanan et al. [18] study Mondrian forests (MFs) as BO surrogates: MFs extend decision forests with a hierarchical Gaussian prior to yield calibrated uncertainty and efficient online updates; on four BO benchmarks (Branin, Hartmann–6, SVM-grid, LDA-grid) MF–UCB matches or exceeds SMAC (RF), and on a large flight-delay regression task MF attains better NLPD than a strong GP approximation (rBCM). The study, however, reports no wall-clock or memory profiles (compute is discussed via O(N log N ) vs. O(N 3 ) arguments rather than measurements), and the BO comparison is against Random Forest (RF) rather than GP. In a representative low-data setting, Tom et al. [19] runs discrete-library BO across six small molecular datasets (150–1,870 candidates) using GP surrogates (including a Tanimoto-kernel variant), NGBoost, SNGP, BNN, and GNNGP, but reports only optimization/calibration results without wall-clock or memory measurements. Likewise, Bayesian model averaging frameworks that integrate non-GP surrogates such as BART and BMARS demonstrate strong optimization outcomes on practical, discrete candidate sets [20], again without compute or memory profiling. An exception on the compute side is the DNGO study of Snoek et al. [21], which includes a direct runtime-scaling plot against a GP baseline (Spearmint) on the Hartmann–6 benchmark—showing roughly linear growth for DNGO versus much steeper increases for the GP—yet this is limited to a single benchmark and does not report memory. Altogether, the field still lacks a systematic, quantitative comparison that elevates both optimization quality and resource cost to first-class metrics, leaving open the central question of whether alternative surrogates can deliver not only accuracy but also genuine efficiency gains over GP-based BO. At the same time, frugality matters: data centers and AI form an energy-hungry duo—DC usage is large and rising [22, 23] and AI’s footprint keeps growing despite efficiency gains, with inference up to ∼60% [24– 26]. Hardware demand and costs are climbing in parallel—driven by AI workloads and

3

acute GPU scarcity [27]—tightening compute budgets in research settings. Even if scientific ML workloads are modest compared to frontier LLMs, adopting a resourceaware mindset (minimizing wall-clock time and memory, and avoiding unnecessary GPUs/HPC) should be standard practice. ML/AI is rapidly being deployed in domains central to net-zero goals (renewable-ready power grids [28], transport and logistics [29], buildings and industrial operations [30], land-use and agriculture [31], and accelerated discovery of materials for batteries [32], renewables [33], Cooling/Dehumidification [34] and CO2 emissions reduction [35]) where it can cut emissions; yet a growing discussion warns that escalating AI compute in net-zero research can offset these gains and even delay net-zero timelines [36, 37]. This tension further motivates resourceaware methods in scientific ML. Historically, GP-based BO has been the default; here we demonstrate that alternative surrogates can replace GPs to deliver comparable or better optimization while reducing compute and memory, enabling more efficient sampling. FruBO framework is a redefinition of how BO should be evaluated. We reframe BO as a two-axis problem (balancing sampling efficiency with computational cost) and provide the first systematic, quantitative study across surrogate models that exposes this trade-off. Within this perspective, efficiency is not measured solely by sample or iteration count, but by the joint behavior of optimization quality and computational scaling (time, memory) as acquisitions grow. We benchmark four surrogates: a GP baseline and three scalable alternatives—Random Forests (RF), Bayesian Adaptive Spline Surfaces (BASS) [38], and NGBoost (NGB) [39]. FruBO offers a practical and transparent way to measure and visualize this balance, establishing reference behavior that future BO research can build on. By standardizing this paired evaluation protocol and releasing it as open-source software, we aim to encourage a culture of computeaware benchmarking in scientific machine learning and to guide the selection and development of surrogate models under realistic resource constraints. We evaluate across two classes of datasets. First, eight well-known benchmark functions which offer controlled landscapes with known optima: Rastrigin, Ackley, Schwefel, Michalewicz, Schaffer 7, Styblinski–Tang, Weierstrass, and Expanded Schaffer F6. Second, 9 domain-relevant problems that differ in representation, smoothness, dimensionality, and noise: (i) a approx. 70,000 COF candidates dataset for maximizing the methane deliverable capacity; (ii) a 400,000-candidate moiré-pattern system of six-layer stacks targeting near-uniform pore-size distributions; (iii) an optimal-control task based on LunarLanderContinuous-v3 reformulated for BO by parameterizing action schedules and using episode return as the black-box objective; (iv) a large-scale molecular screening for therapeutic relevance, looking for the molecule with the lowest docking energy among 400,000 candidates from a large-scale Docking Dataset; (v) ∼50,000 modeled high-pressure tank designs to maximize pressure thresholds; (vi) the 130,000 small molecules QM9 dataset to maximize the HUMO-LUOMO gap value; (vii) an approx. 580,000 candidates dataset of a multi-layer approaximation of a c. elegance model towards minimizing its stiffness; (viii) a neural-network hyperparameter optimization on a scientific task; (ix) and an electromechanical damping towards minimizing structure vibrations across multiple frequencies. Within this framework, we find that non-GP surrogates consistently match or surpass GP-based BO in terms of optimization performance (best-found value and recall@100 as a function of samples), while offering markedly superior computational efficiency (best-found value and recall@100 as a function of compute time). The compute

4

profile differs sharply across models: GP training time grows superlinearly with the number of acquisitions and memory accumulates with each update, whereas RF, BASS, and NGB exhibit nearly linear time scaling with effectively bounded memory usage. In practical terms, BO runs with GP surrogates can terminate prematurely on standard workstations due to memory pressure, even in settings where the objective evaluations themselves are inexpensive. By contrast, the non-GP surrogates remain stable and efficient throughout. These findings demonstrate that the surrogate choice is not a secondary modeling detail but a primary determinant of computational scalability. As acquisition counts increase, the frugality-oriented surrogates deliver equal or better optimization quality at a fraction of the computational cost. Motivated by these results, we take a further step and formalize the selection of surrogates through classifier. Using a small set of easily obtainable dataset characteristics, such as dataset size, dimensionality, fractal dimension, and a lightweight estimate of the target variance, we train a multi-output model that predicts the expected ranking of the four surrogates under either compute-time–limited or sample-acquisition–limited regimes. This approach transforms empirical observations into a practical tool: given a new dataset, users can obtain an informed recommendation for the surrogate model most likely to balance optimization performance and computational frugality. In this way, FruBO not only quantifies the trade-offs inherent in surrogate choice but also operationalizes them, providing an accessible and compute-aware BO strategy for researchers working under constrained data or hardware budgets. In summary, this study combines three complementary elements: a systematic comparison of BO surrogate models across diverse problem classes, a compute-aware evaluation protocol that jointly assesses optimization quality and computational scaling, and a lightweight surrogate-recommendation model that capitalizes on the resulting empirical structure. By linking surrogate benchmarking, evaluation, and selection within a single framework, FruBO provides a practical methodology for deploying Bayesian Optimization in settings where both data and computational resources are constrained.

1 Results 1.1 The premise The premise of this work is that researchers typically operate under fixed experimental or computational budgets while seeking high-performing solutions within large, complex design spaces (for example, materials screening under limited simulation time or laboratory experiments constrained by cost). In our setting, this budget is fixed at N = 1000 new acquisitions (e.g., synthesizing a new material and measuring its performance), which is typically much smaller than the total number of candidates in the design space (N ≪ N ′ , where N ′ denotes the full set of available materials or design options). In addition to optimization outcomes, we explicitly track the computational costs incurred during the BO process. For each surrogate, we report: (i) training time as a function of the number of acquisitions, reflecting how wall-clock cost scales over the course of the loop, and (ii) GPU memory consumption as a function of the number of acquisitions, reflecting whether memory usage remains bounded or grows with dataset size. These two dimensions provide a direct measure of scalability and practical feasibility on standard GPU-equipped workstations, complementing the evaluation of optimization performance. The surrogate models are the standard GP, 5

RF, NGBoost and BASS, and more information for their implementation in our work can be found in the methods section.

1.2 Benchmark functions

Fig. 1: Aggregated BO performance across nine benchmark functions. (a) wall-clock training time vs. acquisitions, and (b) GPU memory vs. acquisitions. Curves are averaged over functions and BO seeds; shaded regions indicate variability.

Figure 1 reports the aggregated compute-performance results (wall-clock time and GPU memory) across all eight benchmark functions, each averaged over 20 BO runs. Shaded regions indicate variability across BO seeds. Because all benchmark functions share identical computational characteristics—namely the same dataset size (up to 1,000 acquired samples) and the same input dimensionality (4D)—this aggregation is meaningful and provides a reliable overview of surrogate-model scaling behavior as the number of BO acquisitions increases. As expected, GP-based BO exhibits the steepest growth in compute cost. The wall-clock training time for GP increases superlinearly with the number of acquired samples, consistent with the well-known O(n3 ) complexity of Gaussian Process regression. The GPU footprint of GP also grows approximately linearly, reflecting the increasing cost of kernel-matrix operations. In contrast, the alternative surrogates, RF, NGB, and BASS, display markedly different scaling. Their wall-clock training time increases only linearly with the number of acquisitions, and their GPU memory usage remains minimal throughout the entire BO trajectory. Importantly, the GPU cost of these models does not increase as more samples are acquired, a consistent and striking result that also appears in the real-world datasets analyzed later. These findings highlight a key practical insight: in settings with larger datasets or higher-dimensional feature spaces (beyond the 4D benchmark functions used here), replacing GP with RF, NGB, or BASS can yield substantial computational savings. Their linear time scaling and negligible GPU footprint make them far more suitable for large-scale or high-throughput BO scenarios, without the prohibitive scaling bottlenecks inherent to GP surrogates. The optimization performance of all surrogates is evaluated with respect to two crucial quantities: compute time and number of acquired samples. Both are important, but their relative weight depends on the nature of the oracle that supplies BO with function evaluations. When the oracle is computational, for instance a lightweight 6

model or a fast equation-based solver, the bottleneck is the BO procedure itself, and compute time becomes the dominant metric. Conversely, when the oracle corresponds to an experiment or an expensive simulation, the true bottleneck is the number of samples, since each evaluation may require hours, days, or even weeks, completely overshadowing the cost of BO. To capture performance across both regimes, we report the evolution of (i) the best objective value found and (ii) the identification of the top–100 best solutions, each plotted as a function of compute time and as a function of the number of acquired samples (four plots in total). Because judging which model performs best is not reliable through visual inspection alone, we quantify performance using the area under the curve (AUC); smaller AUC values indicate better overall performance. In the plots, shaded regions illustrate the variability of AUC across BO runs. Figure 2 summarizes the optimization and compute performance across all benchmark functions.

7

Fig. 2: From top to bottom, the eight benchmark functions considered in this study are shown. From left to right, each row displays: (i) a 3D rendering or PCA-based representation of the design space; (ii) the best objective value found as a function of the number of acquired samples; (iii) the best objective value found as a function of compute time; (iv) the number of top–100 global solutions identified as a function of acquired samples; and (v) the same top–100 identification plotted against compute time. Shaded regions denote AUC variability across BO runs, and all curves compare the 8 four surrogate models (GP, RF, BASS, NGBoost) under identical acquisition settings.

To obtain a single metric describing BO compute-time performance, we multiply the AUC of the “best objective vs. time” curve with the AUC of the “top–100 identification vs. time” curve. This product captures both (i) how rapidly a model approaches the global optimum and (ii) how quickly it identifies a large set of high-quality solutions, all measured in wall-clock time. In an analogous fashion, we estimate BO sampleacquisition performance by multiplying the AUC of the “best objective vs. samples” curve with the AUC of the “top–100 identification vs. samples” curve. All resulting metrics are tabulated in the Supplementary Tables 1–4. In terms of BO compute-time performance, RF is the clear winner across all benchmark functions: it reaches the global optimum fastest and identifies a large fraction of the top–100 solutions within minimal compute time. NGB typically ranks second, while GP is last in most cases. Two exceptions occur: on the Styblinski function, BASS ranks second, and on Schaffer7, GP ranks second rather than last. For BO sample-acquisition performance, NGB dominates, ranking first in six out of eight cases. Again, GP is last or second-to-last in the majority of benchmarks. Taken together, these results (summarized in Table 1) reveal a clear and consistent contrast across all benchmark functions: GP-based BO is substantially more expensive in both compute time and GPU memory, yet this overhead does not translate into superior optimization performance. In fact, the opposite trend emerges. RF, BASS, and NGB achieve faster progress toward the optimum and identify high-quality solutions with far fewer samples, as reflected in their markedly better compute-time and sample-acquisition AUC metrics. These models therefore deliver both stronger optimization performance and significantly lower computational cost. This establishes a decisive expectation for the real-case datasets: surrogate models with linear scaling and negligible GPU requirements—RF, BASS, and NGB—are likely to provide clear advantages when BO is applied to larger, higher-dimensional, or experimentally constrained problems. We now turn to these real-world datasets to assess whether these trends persist in practical settings.

1.3 Real-world datasets The benchmark functions constitute well-behaved testbeds whose response surfaces are dictated by strict analytical formulas, reflecting idealized patterns inspired by physics, biology, or finance. In contrast, real-world optimization problems rarely exhibit such structure: their landscapes are irregular, their sampling is uneven, and the underlying phenomena are far more intricate. To assess Bayesian Optimization under these realistic conditions, we assembled a suite of nine demanding, domain-spanning datasets. These datasets arise either from (i) in-house, expert-developed modeling pipelines constructed specifically for this work, or (ii) carefully curated and processed literature datasets that required substantial domain knowledge, custom preprocessing, and new Python tooling. Together, they span bio-inspired mechanical engineering, biomechanics, electromechanics, drug discovery, materials science, physics, robotics control, and machine learning. In-house, expert-developed datasets • Moiré multilayer materials: geometric patterns generated by superimposed hexagonal layers with controlled rotations. The objective is to maximize the uniformity of the size distribution of the pores created through the rotation of the stacked layers 9

• Bouligand pressure vessel: tens of thousands of composite laminate designs evaluated using a custom finite-element simulation framework for stress-based performance. • Mulit-PZT Semi-active tuned mass damper (SATMD): a structuralvibration dataset targeting multi-modal vibration minimization under semi-active control. • C. elegans worm-indentation mechanics: a mechanics-informed dataset integrating AFM indentation profiles with hyperelastic modeling to quantify drug-induced softening. • Lunar Lander control-sequence dataset: a 50-step discrete action-sequence space defining a black-box control problem through episode-return evaluations. • Neural-network hyperparameter tuning: a wide grid of MLP architectures and training parameters evaluated for predictive performance. Curated and processed literature datasets • Ro4 molecular docking: a curated subset of the ultralarge Ro4 chemical space, processed with RDKit descriptors extraction. • COFs dataset: Monte Carlo simulation data of covalent-organic frameworks for high methane deliverable capacity. • QM9 quantum-chemistry dataset: B3LYP-computed molecular properties for small organic compounds, processed with RDKit descriptors extraction, towards maximizing the LUMO-HOMO gap. Information on the construction, processing, and modeling pipelines for all nine datasets is provided in detail in the Methods section. First of all, it is worth noting that in all real-case datasets the compute-time complexity and GPU-memory scaling as a function of the number of acquired samples follow the same pattern observed for the benchmark functions: Gaussian Processes exhibit a characteristic N 3 growth in compute time and an approximately linear increase in GPU memory. In contrast, RF, NGBoost, and BASS maintain nearly constant GPU memory usage and display linear compute-time scaling throughout the BO trajectory. All corresponding plots—GPU memory and compute time as a function of sample acquisition—are provided in the ESI (Supplementary Figure 3). This again highlights the point raised in the Introduction: GPs impose a substantial computational burden on Bayesian Optimization, making them increasingly strenuous to train as the number of samples grows, whereas the alternative surrogates remain lightweight and scalable.

10

Fig. 3: From top to bottom, the eight real-case data considered in this study are shown. From left to right, each row displays: (i) a 3D rendering or PCA-based representation of the design space; (ii) the best objective value found as a function of the number of acquired samples; (iii) the best objective value found as a function of compute time; (iv) the number of top–100 global solutions identified as a function of acquired samples; and (v) the same top–100 identification plotted against compute time. Shaded regions denote AUC variability across BO runs, and all curves compare the four surrogate models (GP, RF, BASS, NGBoost) under identical acquisition settings.

11

The analysis that follows examines whether these computational observations also translate into improved BO performance, both in terms of compute time and sample efficiency. As before, the results are presented using the same structure as for the benchmark functions: plots of BO compute-time performance and BO sampleacquisition performance for all four surrogate models, with the area under the curve (AUC) serving as a clear and comparable performance metric across datasets. Figure 3 summarizes the optimization and compute performance across all real case datasets. Across the nine real-case datasets, Gaussian Processes exhibit consistently weak performance in both BO compute-time and sample-acquisition metrics. GP ranks first in compute-time performance only twice (the Moire and worm datasets) and achieves the top rank in sample-efficiency only once, again in the worm case, where the landscape happens to align well with GP smoothness assumptions. In all remaining datasets, GP appears as second-to-last or last, highlighting its limited practical competitiveness once computational cost and scalability are taken into account. By contrast, Random Forest and NGBoost emerge as the strongest overall surrogates: RF attains the best compute-time performance in five out of eight datasets, reflecting its linear scaling and robustness, whereas NGBoost dominates sample-efficiency, ranking first in five datasets. BASS performs well in selected cases but is less consistently dominant across domains. All results are summarized in Table 1. The values for all the AUC for all cases and surrogates can be found in All resulting metrics are tabulated in the Supplementary Tables 5–8. Table 1: Ranking of BO surrogate models (1 = best, 4 = worst) in terms of computetime and sample-acquisition performance across benchmark functions and real-case datasets. Compute-time performance Task / Dataset

1

2

3

Sample-acquisition performance

4

1

2

3

4

NGB NGB NGB NGB NGB BASS RF RF

RF BASS RF RF RF NGB BASS NGB

GP RF BASS BASS BASS RF NGB GP

BASS GP GP GP GP GP GP BASS

GP RF BASS RF GP BASS GP RF GP

RF GP GP GP RF RF BASS GP RF

NGB BASS NGB BASS BASS NGB RF NGB BASS

Benchmark Functions expschaffer6 rastrigin Michalewicz Ackley Schwefel Stybilinski Weierstrass Schaffer7

RF RF RF RF RF RF RF RF

NGB NGB NGB NGB NGB BASS NGB GP

BASS BASS BASS BASS BASS GP BASS NGB

GP GP GP GP GP NGB GP BASS

Real-Case Datasets Moire pressure_vessel MOFs Ro4 docking lunar c elegance NN tuning QM9 multi_PZT

GP RF RF BASS RF GP RF BASS RF

BASS NGB BASS RF GP RF NGB RF GP

RF GP NGB NGB BASS NGB GP GP NGB

12

NGB BASS GP GP NGB BASS BASS NGB BASS

BASS NGB RF NGB NGB GP NGB BASS NGB

1.4 Building a recommendation system for surrogate models in BO The left panel of Fig. 4 summarizes the Bayesian Optimization (BO) evaluation pipeline developed and analyzed throughout this work. Across both benchmark functions and real-world datasets, Gaussian Processes (GPs) consistently emerge as the least practical surrogate: their cubic scaling in training time and steadily increasing memory footprint impose a substantial computational burden, while their sampling efficiency rarely compensates for this cost. In contrast, RF, NGBoost, and BASS achieve comparable or superior optimization performance with markedly lower computational overhead. These results highlight a key insight: no surrogate model is universally optimal, and defaulting to GP-based BO is seldom justified across diverse problem settings. Motivated by this observation, we move beyond retrospective comparison and introduce a surrogate-model recommendation strategy. The core idea is to use a small set of dataset characteristics that are inexpensive to obtain, yet informative of surrogate performance. Specifically, we consider four features: dataset size (number of samples), dimensionality (number of input features), fractal dimension of the input space, and the variance of the target property. The first three features are extracted directly from the dataset descriptors, while the target variance is estimated from a short initial BO phase, using the first 200 acquisitions with GP as a surrogate. Based on these features, we train two multi-output classifiers that predict the ranking of the four surrogate models considered in this work (GP, RF, NGBoost, and BASS). One classifier targets compute-time performance, while the other targets sample-acquisition performance. Together, these models act as surrogate recommenders, enabling the selection of an appropriate BO surrogate tailored to the computational or experimental constraints of a given problem. The right panel of Fig. 4 illustrates the complete workflow used to construct these classifiers, transforming empirical BO benchmarking results into a practical, dataset-aware recommendation system. To evaluate the performance of the surrogate-ranking classifier, we adopt a 5fold cross-validation scheme. Each fold is assessed using the Normalized Discounted Cumulative Gain (nDCG) [40], a standard ranking metric that measures the quality of a predicted ordering relative to an ideal ranking, while placing greater emphasis on correctly identifying higher-ranked (more relevant) items. The nDCG at position p is defined as nDCGp =

DCGp , iDCGp

(1)

where DCGp denotes the Discounted Cumulative Gain of the predicted ranking up to position p, and iDCGp is the corresponding gain for the ideal (perfectly ordered) ranking. The Discounted Cumulative Gain is given by DCGp =

p X 2reli − 1 , log2 (i + 1) i=1

where reli is the relevance score of the item at rank i in the predicted ordering.

13

(2)

The ideal Discounted Cumulative Gain is computed analogously by sorting items according to their true relevance: |RELp |

iDCGp =

X i=1

2reli − 1 . log2 (i + 1)

(3)

Here, p denotes the number of items in the ranking, and |RELp | is the number of relevant items considered. For each fold, we report: (i) the mean nDCG across test datasets in that fold, and (ii) the distribution of per-dataset nDCG values aggregated across all folds. This analysis is performed separately for rankings derived from sample-efficiency performance and from computational-efficiency (time-based) performance, allowing us to examine whether surrogate selection can be predicted consistently across different notions of frugality.

Fig. 4: a): our four surrogates evaluation scheme in this work through 8 benchmark functions and 9 real-case datasets; b)Building the classifier that recommends the best surrogate given a dataset

In 5 we show the performance evaluation of our classifier.

14

Fig. 5: Mean nDCG per fold (top) and per-dataset nDCG distributions aggregated across folds (bottom) are shown for rankings derived from sample-efficiency performance (left) and computational-efficiency performance (right).

The results demonstrate that the meta-classifier achieves consistently high nDCG values across folds for both performance criteria. For sample-efficiency rankings (left panels), the mean nDCG per fold remains high, indicating that the predicted surrogate orderings closely match the true rankings derived from BO performance. The per-sample nDCG distribution further shows that most datasets achieve nDCG values near 1, implying that even when the full ranking is not predicted exactly, the classifier typically identifies the most relevant surrogates correctly—particularly at the top of the ranking. A similar, and in fact slightly stronger, trend is observed for computational-efficiency rankings (right panels). Here, the mean nDCG values are consistently high across all folds, and the per-sample distribution is even more concentrated near nDCG ≈ 1. This suggests that surrogate performance with respect to computational cost exhibits more regular, learnable structure across datasets, making it particularly amenable to meta-learning. Overall, these results indicate that the proposed meta-BO framework can reliably infer near-optimal surrogate rankings from simple dataset characteristics, with especially strong performance when targeting computational frugality. The use of nDCG highlights that the classifier is effective not only at exact ranking prediction but, more 15

importantly, at correctly prioritizing the most suitable surrogate models for a given optimization task.

2 Discussion Bayesian Optimization is often presented as a sample-efficient alternative to bruteforce exploration, yet in practice its computational overhead is rarely treated as a first-class concern. In this work, we show that this omission is consequential. Across a diverse set of benchmark functions and domain-relevant datasets, surrogate choice fundamentally determines not only optimization performance but also whether BO remains computationally viable as acquisitions accumulate. Our results demonstrate that the long-standing default Gaussian Processes frequently represents the least practical option once wall-clock time and memory usage are explicitly accounted for. A consistent pattern emerges across all experiments: GP-based BO incurs a rapidly increasing computational and memory burden, while offering little or no compensating advantage in optimization quality. In both benchmark and real-world problems, GP training time grows superlinearly and GPU memory usage accumulates steadily, often becoming the dominant bottleneck even when objective evaluations themselves are inexpensive. By contrast, Random Forests, NGBoost, and BASS exhibit near-linear scaling in compute time and effectively bounded memory usage, allowing BO to proceed smoothly over hundreds to thousands of acquisitions on standard hardware. Crucially, this improved frugality does not come at the expense of optimization success. On the contrary, the non-GP surrogates consistently match or surpass GPbased BO in identifying optimal or near-optimal solutions, both in terms of bestfound values and recovery of top-performing candidates. These findings challenge the widespread assumption that GP smoothness and calibrated uncertainty inherently translate into superior BO performance. Instead, they suggest that, for many realistic design spaces, the computational structure of the surrogate model matters at least as much as its theoretical optimality guarantees. While GP-based BO performs poorly on average, no single alternative surrogate dominates across all datasets and performance regimes. Random Forests excel in compute-time–limited settings, while NGBoost often achieves superior sample efficiency; BASS provides strong performance in selected domains with complex, nonstationary structure. This heterogeneity reinforces a central conclusion of this study: there is no universally optimal surrogate model for Bayesian Optimization. This observation has important practical implications. In experimental science, where sample acquisition is expensive, optimizing for sample efficiency may be paramount. In simulation-driven or data-rich settings, compute time and memory footprint may dominate instead. Treating BO as a one-size-fits-all procedure, by defaulting to GP surrogates, obscures these trade-offs and can lead to unnecessarily costly optimization pipelines. Our results argue for a shift in perspective: surrogate selection should be viewed as a dataset-dependent decision rather than a fixed design choice. Motivated by this insight, we extend FruBO beyond retrospective benchmarking and introduce a data-driven surrogate recommendation framework. Using a small set of inexpensive dataset characteristics (dataset size, dimensionality, fractal dimension, and a lightweight estimate of target variability) we show that it is possible to predict

16

surrogate rankings with high fidelity. The strong nDCG scores achieved by our multioutput classifiers indicate that surrogate performance exhibits learnable structure across datasets, particularly when computational efficiency is the target. This recommendation layer converts empirical benchmarking into a practical decision tool. Instead of relying on costly pilot studies across multiple surrogates, FruBO enables surrogate selection a priori, aligned with either compute-time–limited or sample-limited regimes. In doing so, it moves beyond comparison and embeds computational frugality directly into the BO workflow. The implications extend beyond Bayesian Optimization. As machine learning becomes integral to scientific discovery, computational cost and energy usage can no longer be treated as secondary concerns. Even modest workloads, when repeated across models and optimization loops, accumulate substantial overhead. Our results indicate that compute-aware evaluation protocols can deliver immediate efficiency gains without compromising optimization quality. By reframing BO evaluation along two axes (optimization quality and computational frugality) this work provides a template for future methodological studies in scientific machine learning. We anticipate that similar paired evaluations will be increasingly necessary as the community confronts tighter compute budgets, growing environmental concerns, and the need to democratize advanced optimization tools beyond specialized HPC environments. Several limitations of the present study point to directions for future work. First, while we consider four representative surrogates, the landscape of scalable probabilistic models continues to expand, including sparse GP variants, neural surrogates, and hybrid ensembles. Extending the FruBO framework to include these models would further enrich the recommendation space. Second, our dataset-characteristic features are intentionally simple; incorporating richer descriptors of landscape structure or noise could further improve surrogate selection, particularly in small-data regimes. Finally, while our recommendation system focuses on surrogate choice, future extensions could jointly recommend acquisition functions or dynamically adapt surrogate models during the BO loop. Despite these limitations, the central conclusion is clear: Bayesian Optimization need not be computationally heavy to be effective. By abandoning the default reliance on Gaussian Processes and embracing scalable alternatives guided by data-aware recommendations, BO can become a genuinely frugal and widely applicable tool for scientific discovery.

3 Methods FruBO is implemented on top of the BoTorch framework [41], which provides the core Bayesian Optimization abstractions and Gaussian-process baselines. FruBO extends this foundation by systematically integrating and benchmarking alternative, scalable surrogate models beyond standard GP-based BO, including Random Forests, NGBoost, and Bayesian Adaptive Spline Surfaces (BASS). The framework is designed to be modular, allowing users to easily incorporate additional surrogate models and evaluate them under unified acquisition, evaluation, and logging protocols. The full implementation is publicly available; see the Code Availability section.

17

3.1 Surrogate models Surrogate Models We evaluate four surrogate models within our Bayesian Optimization framework, chosen to balance predictive performance and computational efficiency. Gaussian Processes (SingleTaskGP). We adopt the standard exact GP implementation in BoTorch [41], using automatic relevance determination (ARD) kernels and outcome standardization. GPs remain the canonical choice in BO due to their principled uncertainty quantification, though their training cost scales cubically with the number of observations and memory usage increases linearly. Random Forest (RF) surrogate. Random Forests (RFs) are nonparametric ensemble models that combine multiple regression trees to approximate complex, nonlinear response surfaces. We include RF as a frugal and well-established baseline surrogate, widely adopted in modelbased optimization and algorithm configuration [42]. RFs naturally handle nonlinear, nonstationary, and mixed-type inputs without requiring kernel definitions or feature scaling, and their training cost scales nearly linearly with the number of samples and trees, making them particularly suitable for iterative Bayesian Optimization loops. Compared to Gaussian Processes, RFs avoid the O(n3 ) covariance inversion bottleneck while providing variance estimates from the ensemble dispersion, albeit with coarser calibration. In contrast to more sophisticated surrogates such as BASS and NGBoost, RFs are computationally lighter, require minimal hyperparameter tuning, and are robust to outliers and discontinuities. We use the standard RandomForestRegressor implementation from scikit-learn with default parameters. Uncertainty is estimated from the variance across individual trees. While RF is less theoretically grounded in its uncertainty quantification compared to GPs, it remains computationally efficient, robust to noise, and scales well to large datasets. Their interpretability, mature implementations, and strong empirical record in sequential model-based optimization (SMBO) frameworks (e.g., SMAC [42]), makes RF a good alternative to GP as scalable BO surrogates [43]. Bayesian Adaptive Spline Surfaces (BASS). BASS [38] is a nonparametric Bayesian regression model based on adaptive spline bases with posterior inference via reversible-jump MCMC sampling. It extends the Bayesian Multivariate Adaptive Regression Splines (BMARS) framework of Denison et al. [44] by supporting both continuous and categorical inputs, scalar or functional outputs, and by incorporating more flexible priors and Reversible Jump Markov Chain Monte Carlo (RJMCMC) strategies for exploring model space. BASS is also equipped with analytical Sobol decompositions, making it directly useful for sensitivity analysis. These enhancements make BASS more versatile and scalable than BMARS for high-dimensional or mixed-type design spaces. We include BASS among our surrogate models because it provides a computationally frugal, nonstationary-capable, and uncertainty-aware alternative to Gaussian Processes. Its adaptive spline basis allows the model to vary smoothness across the input space without assuming stationarity, capturing heterogeneous regions that often 18

appear in physical and benchmark functions. The RJMCMC inference mechanism automatically adjusts model complexity by adding or pruning basis functions, while shrinkage priors control overfitting as the BO loop progresses. Compared to BMARS, BASS offers improved scalability, parallelization options, and more efficient posterior exploration through parallel tempering. It natively handles categorical predictors and functional responses, and its full Bayesian formulation yields structured uncertainty estimates readily usable by acquisition functions. Together, these features make BASS particularly suited to frugality-oriented Bayesian Optimization, where flexibility, interpretability, and efficient uncertainty quantification are required under limited evaluation budgets. Across domains, BASS has demonstrated state-of-the-art predictive accuracy and efficiency. In building-energy modeling, multi-output BASS achieves CV(RMSE) < 0.005 with R2 ≈ 1, outperforming deep-learning baselines while preserving correlations among multiple outputs, and its analytical Sobol’ implementation reproduces sampling-based rankings with 46×–148× faster sensitivity computation [45, 46]. In materials and nuclear-forensics applications, BASS enables accurate inverse prediction of plutonium processing conditions—on average within one standard deviation of the true experimental settings [47]—and shows top performance in predicting detonation metrics such as detonation velocity (Vdet ) and detonation pressure (Pdet ) among several machine-learning models [48]. These results confirm BASS as a robust and high-performing surrogate for complex, high-dimensional, and mixed-type problems. However, to our knowledge, it has not yet been implemented or evaluated within a full Bayesian Optimization framework, which motivates its integration and benchmarking in this work. BASS was originally developed as an R package [38] (see manual at https://cran. r-project.org/web/packages/BASS/BASS.pdf) but has recently been ported to Python as pyBASS (https://github.com/lanl/pyBASS/tree/main), which we employ in this work. We follow settings close to the defaults, but with an even more frugal choice of burn-in: whereas the default is nmcmc = 10000 and nburn = 9000 [38], we use nmcmc = 10000 and nburn = 9900. This reduces the number of effective posterior samples, but substantially lowers computational cost when used inside a repetitive BO pipeline. In this sense, our configuration prioritizes frugality over extracting the very highest predictive performance. nmcmc nburn thin w1, w2 maxInt maxBasis g1, g2 s2_lower h1, h2 a_tau b_tau verbose

= 10000 % total MCMC iterations = 9900 % burn-in iterations (default 9000) = 1 = 5.0, 5.0 = 2 = 1000 = 1.0, 1.0 = 0.0 = 10.0, 10.0 = 0.5 = 1.0 = True

19

Natural Gradient Boosting (NGBoost). NGBoost is a modular probabilistic boosting framework that uses the natural gradient to learn full predictive distributions with stable training dynamics and strong scalability, offering both accurate point predictions and well-calibrated uncertainty estimates [39]. Empirically, in low-data chemical design and Bayesian Optimization (BO) benchmarks, NGBoost achieved among the top fractions of hits (e.g., Delaney dataset 0.959, Freesolv 0.953) and exhibited strong early-stage optimization and calibration performance, often matching Gaussian Processes around ∼100 samples [19]. In engineering optimization, a BO–NGBoost surrogate for tunnel deformation prediction attained a test R2 ≈ 0.92 and outperformed Random Forest, XGBoost, LightGBM, GRU, and LSTM models, enabling multi-objective improvements up to ∼56–62% [49]. Taken together, these results confirm NGBoost as a robust and high-performing surrogate capable of providing reliable probabilistic predictions under limited data, motivating its inclusion and systematic benchmarking within our Bayesian Optimization framework. In principle, NGBoost defaults to a relatively conservative configuration: n_estimators=500, learning_rate=0.01, minibatch_frac=1.0, col_sample=1.0, verbose=True, verbose_eval=100, tol=1e-4, random_state=None, validation_fraction=0.1

However, in practice, Duan et al. (2020) also demonstrate that a larger learning rate (0.1) is appropriate for large datasets. We therefore follow this latter choice and at the same time reduce the number of estimators to 100. This configuration is more computationally economical, making it suitable for iterative use within BO, where many surrogate fits are required. While this choice may sacrifice some predictive accuracy compared to the full defaults, it aligns with our emphasis on frugality and adaptability to medium-scale compute resources. n_estimators = 100 % reduced from default 500 learning_rate = 0.1 % increased from default 0.01 random_state = 42 natural_gradient = True score = "LogScore" verbose = False early_stopping_rounds = None sigma_floor = 1e-2 sigma_cap = 1e2 Together, these four surrogates allow us to benchmark classical GP-based BO against alternative ensemble and boosting approaches that bypass the cubic scaling bottleneck of GPs while offering competitive or superior sample efficiency. Importantly, for both BASS and NGBoost we deliberately adopt computationally frugal settings, as our goal is not to maximize surrogate accuracy in isolation, but to evaluate their effectiveness within the demanding, repetitive context of Bayesian Optimization.

20

3.2 Datasets 3.2.1 COFs (methane uptake and deliverable capacity). Dataset size: 69,839; dimensionality: 20; target: deliverable capacity; objective: maximize. Mercado et al. [50] reported ∼70,000 COFs with Monte Carlo simulations for CH4 uptake and deliverable capacity. Here we target deliverable capacity as separate objectives within our BO framework.

3.2.2 QM9 Dataset and Molecular Descriptor Extraction. Dataset size: 133,885; dimensionality: 46; target: HOMO-LUMO gap; objective: maximize. We employed the QM9 dataset introduced by Ramakrishnan et al. [51], which provides B3LYP/6-31G(2df,p)-level quantum-chemical properties for 133,885 stable organic molecules composed of C, H, N, O and F. The dataset is distributed as extended .xyz files, where each entry contains not only the atomic coordinates but also all computed scalar properties (including orbital energies, polarizabilities, heat capacities, thermochemical quantities, and the HOMO–LUMO gap). Crucially for our workflow, each file additionally stores two SMILES representations (the original GDB17 SMILES and a geometry-relaxed SMILES), enabling deterministic reconstruction of each molecule within RDKit for descriptor calculation. We did not use the full RDKit descriptor library but instead retained a subset of 43 physico-chemical descriptors appropriate for QM9. These include global size and composition features (e.g. MolWt, ExactMolWt, NumValenceElectrons), topological and complexity indices (Chin, Kappan, BalabanJ, Ipc), polarity-related descriptors (TPSA, MolLogP, EState indices), and ring- and functionality-count descriptors. This reduced set avoids redundancy, excludes descriptors that are identically zero within the restricted chemical space of QM9, and captures the major structural factors that influence frontier-orbital energetics. For the present analysis, the target property was the HOMO–LUMO gap (gap) provided directly in the QM9 files.

3.2.3 Ro4 docking dataset. Dataset size: 397,104; dimensionality: 210; target: binding energy; objective: minimize. We further incorporated a large-scale dataset originating from the recent work of Lüttens et al. [52], who combined molecular docking with machine learning to accelerate ultralarge virtual screening campaigns. The dataset, deposited on Zenodo (record 10.5281/zenodo.7953917), contains docking scores for a multi-billion “rule-offour” (Ro4) chemical subspace derived from the Enamine REAL library. Ro4 molecules satisfy the property constraints of molecular weight <400 Da and cLogP <4, and were docked against a set of therapeutically relevant protein targets. The primary target considered here is the A2A adenosine receptor (A2AR), for which the dataset provides per-compound docking energies (kcal/mol). We adopt these docking scores as the optimization objective, treating them as a black-box property to be maximized in the BO loop. Because the raw dataset comprises billions of molecules, we reduced its size for tractability. Specifically, we first sampled approximately 5% of the available entries, and then retained 60% of this subset after filtering out compounds with docking energies

21

worse than 5000 kcal/mol. This yielded a final working dataset of 397,104 molecules. Each compound was represented by the full set of 208 RDKit molecular descriptors (physicochemical, topological, and electronic features), resulting in a 209-dimensional dataset once the docking score was included as the target property. This preprocessing pipeline produced a compact yet chemically diverse benchmark that preserves the challenging distribution of docking scores while remaining computationally manageable for surrogate modeling within FruBO.

3.2.4 Lunar Lander Benchmark dataset Dataset size: 500,000; dimensionality: 20; target: episode return; objective: maximize. Dataset motivation and scope To evaluate our sampling and optimization strategies on a real control problem, we implemented a discrete action–sequence benchmark based on the LunarLander-v3 environment from Gymnasium (OpenAI Gym). Our formulation is inspired by the active optimization setting introduced Wei et al. [53], where the classical reinforcementlearning control task is reformulated as a black-box design-space optimization problem. Instead of learning a policy, the objective is to discover the best open-loop sequence of discrete actions that maximizes the total episode reward. Design space The environment exposes a discrete action space A = {0, 1, 2, 3} = {no-op, left engine, main engine, right engine}, and each candidate solution is represented as a fixed-length action vector. While Wei et al. employ a horizon of 100 steps (resulting in a 100-dimensional design space with 4100 possible action sequences), we use a compressed horizon of 50 steps (K = 50), which preserves the structure of the optimization problem while reducing temporal complexity and computational cost. Each candidate action sequence a = (a1 , a2 , . . . , a50 ),

at ∈ A,

is executed deterministically from a fixed initial state using a controlled random seed. During execution, the environment applies its native reward shaping (distance-totarget minimization, leg-contact bonuses, fuel usage and crash penalties), and the cumulative sum of these rewards is recorded as a scalar episode_return. For dataset construction, we uniformly sample sequences from A50 and evaluate them in parallel, producing rows of the form: (a1 , a2 , . . . , a50 , episode_return). This results in a discrete, high-dimensional, black-box design space analogous to the formulation of Wei et al., but with a reduced horizon (50 vs. 100 actions) enabling more efficient benchmarking of frugal optimization methods.

22

Objective for BO: Reward function. In the LunarLander-v3 environment, the objective of each action sequence is to maximize the cumulative return obtained from a shaped reward function that guides the lander toward a soft, upright touchdown. Let the lander state at time t be (xt , yt , ẋt , ẏt , θt , θ̇t , c1,t , c2,t ), where (xt , yt ) is the position relative to the landing pad, (ẋt , ẏt ) are translational velocities, θt and θ̇t the orientation and angular velocity, and c1,t , c2,t ∈ {0, 1} indicate left- and right-leg ground contact. The instantaneous reward used by LunarLander-v3 is q q rt = − 100 x2t + yt2 − 100 ẋ2t + ẏt2 − 100 |θt | − 10 |θ̇t | (4) + 10 c1,t + 10 c2,t + Rengine (at ), where the engine-use penalty is   −0.30,   Rengine (at ) = −0.03,    0,

at = main engine, at ∈ {left engine, right engine}, at = no-op.

A terminal bonus or penalty is applied upon episode termination: (

rterminal =

+100,

successful landing with both legs down,

−100,

crash or unstable contact.

The overall objective returned for an action sequence a = (a1 , . . . , aK ) is therefore episode_return(a) =

T X

rt + rterminal ,

(5)

t=1

where T ≤ K is the first time step at which the episode terminates (success, crash, or timeout). This scalar quantity is used as the target value in our Bayesian Optimization experiments. Figure 6 schematically summarizes the reward construction used in the LunarLander benchmark. Starting from a fixed initial state, a candidate open-loop action sequence induces a trajectory of lander states. At each time step, the instantaneous reward rt is formed by combining state-dependent shaping terms—penalizing distance to the landing pad, velocity, tilt, and angular rate—with action-dependent engine-use penalties. Upon landing, crash, or timeout, a terminal reward rterminal is applied. The final objective value, episode_return, is obtained by summing all step rewards and the terminal contribution, as defined in Eqs. (4)–(5).

23

Fig. 6: Schematic illustration of the LunarLander reward construction. An open-loop action sequence generates a trajectory from an initial to a terminal state. At each time step, state-dependent shaping terms and action-dependent engine-use penalties contribute to the step reward rt , while a terminal reward is added upon landing, crash, or timeout to yield the final episode return.

3.2.5 MLP Hyperparameter Grid Dataset Dataset size: 8,640; dimensionality: 18; target: R2 ; objective: maximize. To construct a supervised dataset for neural-network hyperparameter tuning, we systematically evaluated a feed-forward multilayer perceptron (MLP) on the Lunar_landing regression task. We defined a discrete grid over architectural and training hyperparameters as follows: hidden-layer architectures with depth ∈ {2, 3, 4, 5, 6} and uniform width per layer ∈ {8, 16, 64, 128} (20 architectures in total); activation function ∈ {relu, tanh, logistic}; solver fixed to adam; L2 regularization strength α ∈ {10−6 , 10−5 , 10−4 , 10−3 , 10−2 , 10−1 }; mini-batch size ∈ {32, 64, 128, 256}; and initial learning rate ∈ {10−6 , 10−5 , 10−4 , 10−3 , 10−2 , 10−1 }. Early-stopping controls were fixed to max_iter=50, tol=5×10−2 , and n_iter_no_change=10. This grid results in 8640 unique hyperparameter configurations, each of which was trained using standardized inputs (StandardScaler) and evaluated using the coefficient of determination (R2 ), computed via two ShuffleSplit folds with 20% training and 80% validation data under a fixed random seed. The grid evaluation was executed in parallel using joblib with 15 worker processes and one BLAS thread per process to avoid oversubscription. All evaluated configurations, together with their validation R2 and training time, were stored in a comprehensive CSV file, from which we derived a compact “slim” dataset (nn_grid_lunar_score.csv) containing only the hyperparameters (including

24

numeric summaries of architecture depth and width, and one-hot encodings of activation and solver choices) and the corresponding validation score, with duplicate settings removed.

3.2.6 Aldicarb-Induced Softening: Mechanics-Informed Inverse-Design Dataset Dataset size: 579,681; dimensionality: 3; target: mse; objective: minimize. Dataset motivation and scope Neuromuscular agents such as aldicarb have long served as model compounds in C. elegans pharmacological assays, traditionally probed through behavioral readouts such as motility or paralysis [54]. Recent empirical and in-silico work, however, has revealed that aldicarb produces a pronounced mechanical phenotype: instead of hypercontraction, the dominant effect is a bulk softening of the organism arising primarily from depressurization of the pseudocoelom [55]. This process leads to a stiffness reduction of approximately 66% relative to BDM-treated controls and correlates with volumetric shrinkage and decreased internal hydrostatic pressure. As a result, aldicarb acts as a multi-tissue modulator of organismal biomechanics rather than a purely neuromuscular agent. Building upon this characterization, we constructed a mechanicsinformed dataset representing the limiting mechanical regime of aldicarb-induced softening. The dataset is designed specifically for benchmarking surrogate-modelbased Bayesian Optimization methods on a large-scale, nonlinear inverse problem. The underlying experimental data consist of atomic force microscopy (AFM) force– displacement (F –δ) curves from worms exposed to aldicarb in the high-softening regime. These measurements provide the ground-truth mechanical response against which candidate material parameter sets are evaluated. Following the multilayered anatomical structure of C. elegans, three effective Young’s moduli are considered: E1 (cuticle),

E2 (muscle layer),

E3 (pseudocoelomic region).

These parameters define the elastic state of the organism in a simplified, axisymmetric indentation model, consistent with Neo-Hookean hyperelasticity, incompressibility (ν = 0.5), and previously validated finite-element representations [56]. Figure 7 illustrates the three-layer mechanical model of C. elegans used to relate tissue-level elastic properties to AFM indentation force–displacement responses.

25

Fig. 7: Schematic cross-section of the three-layer indentation model of C. elegans used in this work. An AFM indenter applies a normal force F and induces a displacement δ on the worm body. The organism is modeled as a multilayered structure comprising an outer cuticle layer (E1 ), an intermediate muscle layer (E2 ), and an inner pseudocoelomic region (E3 ). These effective Young’s moduli parameterize the tissue-scale mechanical response and define the inverse-design space explored via Bayesian Optimization.

Design space To build a comprehensive benchmark for BO, we enumerated a full Cartesian grid in the space of layer-wise elasticities using experimentally informed bounds: E1 ∈ [30, 250],

E2 ∈ [300, 3000],

E3 ∈ [30, 500].

A uniform discretization step of ∆E = 5 kPa was applied to all axes. After filtering nonphysical combinations, the resulting design space contained 579,681 unique elasticity triplets. Each triplet (E1 , E2 , E3 ) represents a candidate mechanical phenotype of the worm. For each configuration, we use a nonlinear multilayer stiffness model as presented by Rekatsinas et al. [56]. This model computes the predicted indentation force profile Fpred (δ) over the experimentally measured displacements. The model combines: 1. Hertzian and post-Hertzian indentation regimes with displacement-dependent power-law exponents; 2. beam-theoretic bending stiffnesses of the cuticle, muscle, and gut layers through geometry-specific moments of inertia; 3. hyperelastic correction factors capturing nonlinear deformation at large indentation depths; 4. a series coupling between indenter stiffness and tissue stiffness.

26

This generates a deterministic and physically well-grounded forward map (E1 , E2 , E3 ) 7−→ Fpred (δ). Objective for BO: Force-Matching Error. To quantify how well each candidate parameter set reproduces the experimental indentation curve, we compute the mean squared error: N

MSE(E1 , E2 , E3 ) =

2 1 X Fpred (δi ; E1 , E2 , E3 ) − Ftrue (δi ) . N i=1

Each row of the dataset thus contains three input features (E1 , E2 , E3 ) and a single scalar regression target (MSE), forming a high-dimensional, non-convex, anisotropic inverse-design landscape ideally suited for testing BO surrogates and acquisition strategies.

3.2.7 High-pressure vessel Dataset size: 52,272; dimensionality: 5; target: scalarized objective J; objective: minimize. Dataset motivation and scope Although Bouligand (helicoidal) laminates are widely recognized for their ability to twist cracks, diffuse damage, and suppress delamination, no consolidated, machinelearning-ready dataset exists for pressure-retaining CFRP shells. Existing studies highlight the promise of helicoidal architectures—including graded, hemi-symmetric, and discontinuous variants—for improving damage tolerance [5], yet raw data are sparse or unavailable. To support surrogate-model comparison under realistic mechanical conditions, we constructed a large, structured design–response dataset tailored to Bayesian Optimization (BO) for architected composite materials. Design space The Bouligand generator spans five discrete, physically interpretable parameters: starting angle (0:5:175◦ ), inter-ply pitch (5:5:55◦ ), ply count (8:1:40), symmetry flag {0, 1}, and ply-thickness mode {1, 2}. The full enumeration (36 × 11 × 33 × 2 × 2) yields 52,272 unique laminate layups. Figure 8 schematically illustrates the Bouligand composite pressure-vessel architecture and highlights the key design parameters governing the helicoidal laminate stacking sequence, including ply orientation, inter-ply pitch, thickness, and symmetry.

27

Fig. 8: Schematic illustration of the bio-inspired Bouligand composite pressurevessel architecture. The laminate is constructed from multiple plies with gradually rotating fiber orientations, parameterized by a starting angle α, inter-ply pitch p, ply enumeration through the thickness, ply thickness t, and a symmetry flag controlling mid-plane mirroring. Together, these variables define the helicoidal stacking sequence used to generate the pressure-vessel design space evaluated in this work.

Simulation setup For each configuration, a CFRP pressure vessel (cylindrical body with spherical end caps; 165 mm diameter, 350 mm half-length) is simulated under an internal pressure of 62.5 MPa. Abaqus/CAE generates the geometry, symmetry boundary conditions, and S8R/STRI65 meshing. A unidirectional IM7/8551 lamina is modeled as an orthotropic ply with Hashin damage initiation [57]. Field outputs include section-point stresses S11 and S22 , displacements, reaction forces, and failure flags. Objective for BO. To compare surrogate models, we scalarize the multi-objective problem (minimizing S11 , S22 , and thickness) into: min J = 0.48

S11 S22 + 0.40 + 0.12 Thick, Xt Yc

(6)

where (Xt , Yc ) = (2,560 MPa, 185 MPa) follow a maximum-strength criterion [58]. Weighted-sum scalarization [59] is widely used to generate supported Pareto-optimal solutions; here it provides a unified objective for BO-based surrogate comparison.

3.2.8 Moiré Dataset Dataset size: 400,000; dimensionality: 6; target: Jain’s index; objective: maximize. Dataset motivation and scope Moiré patterns arise when two or more periodic meshes, grids, or lattices are superimposed with a slight mismatch in orientation, 28

spacing, or alignment [60]. Figure 9 shows the result of rotation of one layer over another, for six (6) rotation angles. Even if each individual lattice is perfectly regular, their overlaid configuration generates a new, larger-scale interference pattern composed of slowly varying regions of local alignment and misalignment [61]. Small angular offsets or spacing differences accumulate over distance, giving rise to extended “supercells” whose characteristic length scale can be orders of magnitude larger than the primitive lattice [62]. In modern materials research, Moiré patterns serve both as a sensitive probe and a tunable design parameter [63]. Their effects span natural and engineered systems. In van der Waals heterostructures—such as rotated bilayer graphene, transition-metal dichalcogenides, and related 2D materials—Moiré superlattices can radically reshape the electronic band structure, enabling phenomena central to twistronics and quantum devices [64]. Beyond electronics, Moiré geometry can tailor transport and mechanical response in porous membranes [65], photonic structures [66], and magnetic materials [67].

Fig. 9: Moire.

Design space For each independent sample in our dataset, a multilayer hexagonal mesh was constructed using the hexalattice library (https://github.com/alexkaz2/hexalattice). In our study, we employed six layers. The first layer was fixed at a rotation of 0◦ , while the remaining five layers were assigned random rotation angles drawn uniformly from 0◦ to 60◦ . All hexagons were rendered without fill, with a minimum hexagon diameter of one unit. A centrally located window defined the fixed field of view. After rendering, each multilayer structure was rasterized to an RGBA image buffer, cropped, converted to grayscale, and thresholded using Otsu’s global method to isolate mesh linework. Connected components were then identified via standard region-labeling procedures. Only components fully contained within the cropped frame were retained, excluding any touching the image boundary. For each sample, the areas of all valid connected components (i.e., pores) were computed, and a histogram of pore areas was constructed using linearly spaced bins from zero to the maximum pore area observed in that sample. Using this procedure, we generated a total of 400,000 samples, each corresponding to a unique multilayer configuration and its resulting Moiré pore-area distribution. This large-scale dataset provides a high-dimensional, highly nonlinear, and geometrically complex landscape suitable for benchmarking surrogate-based Bayesian Optimization.

29

Objective for BO. For every pore-area distribution, we compute Jain’s fairness index, a classical measure of distributional uniformity: P 2 ( i ai ) J = P 2, n i ai where ai are the pore areas and n is the number of pores. Values close to 1 indicate highly uniform pore sizes, while lower values denote strong heterogeneity. This scalar serves as the BO target, quantifying how rotation-induced Moiré interference influences geometric regularity.

3.2.9 Multi-PZT Semi-Active Tuned Mass Damper (SATMD) Dataset size: 200,000; dimensionality: 7; target: multi-modal vibration; objective: minimize. Semi-active tuned mass dampers (SATMDs) are widely used for vibration mitigation in large, lightweight, and flexible structures [68–70]. An SATMD combines an auxiliary mass with a piezoelectric element connected to a shunted RL network, enabling real-time modification of the system’s effective stiffness and damping properties. By appropriately tuning the resistance and inductance values of the shunt circuit, the device can suppress vibration amplitudes within targeted frequency bands. Determining the optimal tuning parameters, however, is non-trivial, as the SATMD dynamics depend simultaneously on the modal characteristics of the host structure and on the properties of the external excitation. In multi-PZT SATMD configurations, multiple piezoelectric elements are incorporated into the same device, with each element equipped with its own independent shunt circuit. This arrangement introduces multiple anti-resonances into the host system, enhances multi-modal vibration attenuation, and enables effective vibration suppression over a broader frequency range. At the same time, the dimensionality of the tuning problem increases substantially, as multiple resistance and inductance parameters must be optimized concurrently, rendering exhaustive search strategies impractical. To construct the dataset used in this work, a multi-PZT SATMD configuration was evaluated on a representative test structure exhibiting four dominant structural modes at approximately 36 Hz, 47 Hz, 65 Hz, and 159 Hz, as identified from the frequency response function (FRF) of the measured vertical acceleration at the fuselage tip. The SATMD comprised three piezoelectric elements, each connected to its own RL shunt circuit, while an auxiliary mass contributed a fourth anti-resonance. The resulting optimization variables are the auxiliary mass ratio md —defined as the ratio between the damper mass and a reference mass of the host structure—and the six shunt-circuit parameters (R1 , R2 , R3 , L1 , L2 , L3 ) associated with the three piezoelectric elements, yielding a seven-dimensional design space. Objective for Bayesian Optimization. The objective of the optimization is to suppress the system response within frequency regions where the vibration magnitude lies at least 20 dB below the modal peaks. Figure 10 illustrates this objective by comparing the baseline frequency response of the host structure, characterized by four dominant structural modes, with the attenuated response obtained using a multi-PZT SATMD for a representative configuration of the auxiliary mass ratio and shunt-circuit parameters. 30

Fig. 10: Photograph of the experimental test rig used to evaluate the multi-PZT SATMD mounted on a flexible host structure (center), with the damper location highlighted. A schematic representation of the multi-PZT SATMD configuration is shown (inset), illustrating the auxiliary mass and multiple piezoelectric elements connected in series, each equipped with an independent shunt circuit. Frequency response functions of the vertical acceleration at the measurement location are shown for the baseline structure (solid black) and for a representative SATMD configuration (dashed red), demonstrating simultaneous attenuation of the four dominant structural modes within the frequency bands targeted by the BO objective.

Since lower-frequency modes typically induce larger vibration amplitudes, frequencydependent weighting factors were introduced to emphasize mitigation in the most critical bands. The optimal SATMD parameters are obtained by minimizing the following weighted objective function: min f (SATMDmd ,R1 ,L1 ,R2 ,L2 ,R3 ,L3 ) = w1

Z 37 Hz

|ütip,SATMD (f )| df

35 Hz Z 48 Hz

+ w2 + w3

46 Hz Z 68 Hz

|ütip,SATMD (f )| df (7) |ütip,SATMD (f )| df

62 Hz

+ w4

Z 162 Hz

|ütip,SATMD (f )| df.

156 Hz

A total of 2 × 105 distinct SATMD configurations were generated within the prescribed parameter bounds. For each configuration, vibration attenuation was evaluated in the vicinity of the four structural modes, and the weighted objective value was computed accordingly.

31

Acknowledgements. This work was supported by the European Union’s Horizon Europe research and innovation programme under grant agreement No. 101135927 (NOUS). Supporting Information. Data availability.

Supporting Information accompanies this publication.

Datasets are deposited at 10.5281/zenodo.18194022.

Code availability. The BO and evaluation metrics estimation codes underlying this work are freely available for general use under Apache 2.9 LICENSE (https: //www.apache.org/licenses/LICENSE-2.0) and are deposited at https://github.com/ insane-group/FruBO/. Declarations.

The authors declare no competing interests.

Editorial Policies for: Springer journals and proceedings: https://www.springer.com/gp/editorial-policies Nature Portfolio journals: https://www.nature.com/nature-research/editorial-policies Scientific Reports: https://www.nature.com/srep/journal-policies/editorial-policies BMC journals: https://www.biomedcentral.com/getpublished/editorial-policies

Appendix A

Section title of first appendix

An appendix contains supplementary information that is not an essential part of the text itself but which may be helpful in providing a more comprehensive understanding of the research problem or it is information that is too cumbersome to be included in the body of the paper.

References [1] Settles, B.: Active Learning Literature Survey. Technical Report January, University of Wisconsin–Madison, Madison (2009) [2] Shahriari, B., Swersky, K., Wang, Z., Adams, R.P., De Freitas, N.: Taking the human out of the loop: A review of Bayesian optimization. Proceedings of the IEEE 104(1), 148–175 (2016) https://doi.org/10.1109/JPROC.2015.2494218 [3] Chitturi, S.R., Ramdas, A., Wu, Y., Rohr, B., Ermon, S., Dionne, J., Jornada, F.H.d., Dunne, M., Tassone, C., Neiswanger, W., Ratner, D.: Targeted materials discovery using Bayesian algorithm execution. npj Computational Materials 10(1), 1–12 (2024) https://doi.org/10.1038/s41524-024-01326-2 arXiv:2312.16078 [4] Wu, Y., Walsh, A., Ganose, A.M.: Race to the bottom: Bayesian optimisation for chemical problems. Digital Discovery 3(6), 1086–1100 (2024) https://doi.org/10. 1039/d3dd00234a

32

[5] Loutas, T., Oikonomou, A., Rekatsinas, C.: Bio-inspired discontinuous composite materials with a machine learning optimized architecture. Composite Structures 351, 118597 (2025) https://doi.org/10.1016/j.compstruct.2024.118597 [6] Krokidas, P., Gkatsis, V., Theocharis, J., Giannakopoulos, G.: Navigating materials design spaces with efficient Bayesian optimization: a case study in functionalized nanoporous materials. Digital Discovery 4, 3753–3763 (2025) https://doi.org/10.1039/d5dd00237k [7] Rasmussen, C.E., Williams, C.K.I.: Gaussian Processes for Machine Learning. MIT Press, Cambridge, MA (2006) [8] Siemenn, A.E., Buonassisi, T.: Decreasing the Computing Time of Bayesian Optimization Using Generalizable Memory Pruning. In: 2023 IEEE High Performance Extreme Computing Conference, HPEC 2023 (2023). https://doi.org/10.1109/ HPEC58863.2023.10363486 [9] Tani, L., Rand, D., Veelken, C., Kadastik, M.: Evolutionary algorithms for hyperparameter optimization in machine learning for application in high energy physics. European Physical Journal C 81(2), 1–9 (2021) https://doi.org/10.1140/epjc/ s10052-021-08950-y arXiv:2011.04434 [10] Tani, L., Veelken, C.: Comparison of Bayesian and particle swarm algorithms for hyperparameter optimisation in machine learning applications in high energy physics. Computer Physics Communications 294, 108955 (2024) https://doi.org/ 10.1016/j.cpc.2023.108955 arXiv:2201.06809 [11] Santoni, M.L., Raponi, E., De Leone, R., Doerr, C.: Comparison of HighDimensional Bayesian Optimization Algorithms on BBOB. ACM Transactions on Evolutionary Learning and Optimization 4(3) (2024) https://doi.org/10.1145/ 3670683 arXiv:2303.00890 [12] Siemenn, A.E., Ren, Z., Li, Q., Buonassisi, T.: Fast Bayesian optimization of Needle-in-a-Haystack problems using zooming memory-based initialization (ZoMBI). npj Computational Materials 9(1) (2023) https://doi.org/10.1038/ s41524-023-01048-x arXiv:2208.13771 [13] Wei, Y., Zhuang, V., Soedarmadji, S., Sui, Y.: Scalable Bayesian Optimization via Focalized Sparse Gaussian Processes. Advances in Neural Information Processing Systems 37(NeurIPS), 1–25 (2024) [14] Chang, P.E., Verma, P., John, S.T., Solin, A., Khan, M.E.: Memory-Based Dual Gaussian Processes for Sequential Learning. Proceedings of Machine Learning Research 202, 4035–4054 (2023) arXiv:2306.03566 [15] Lin, Y., Gao, B., Tang, J., Zhang, Q., Qian, H., Wu, H.: Deep Bayesian active learning using in-memory computing hardware. Nature Computational Science 5(January) (2024) https://doi.org/10.1038/s43588-024-00744-y [16] Makrygiorgos, G., Ip, J.H.S., Mesbah, A.: Towards Scalable Bayesian Optimization

33

via Gradient-Informed Bayesian Neural Networks. IFAC-PapersOnLine 59(6), 157–162 (2025) https://doi.org/10.1016/J.IFACOL.2025.07.138 arXiv:2504.10076 [17] Li, Y.L., Rudner, T.G.J., Wilson, A.G.: a Study of Bayesian Neural Network Surrogates for Bayesian Optimization. In: 12th International Conference on Learning Representations, ICLR 2024 (2024) [18] Lakshminarayanan, B., Roy, D.M., Teh, Y.W.: Mondrian forests for large-scale regression when uncertainty matters. Proceedings of the 19th International Conference on Artificial Intelligence and Statistics, AISTATS 2016 51, 1478–1487 (2016) arXiv:1506.03805 [19] Tom, G., Hickman, R., Zinzuwadia, A., Mohajeri, A., Sanchez-Lengeling, B., Aspuru-Guzik, A.: Calibration and generalizability of probabilistic models on low-data chemical datasets with DIONYSUS. Digital Discovery 2(3), 759–774 (2023) https://doi.org/10.1039/d2dd00146b arXiv:2212.01574 [20] Lei, B., Kirk, T.Q., Bhattacharya, A., Pati, D., Qian, X., Arroyave, R., Mallick, B.K.: Bayesian optimization with adaptive surrogate models for automated experimental design. npj Computational Materials 7(1), 1–12 (2021) https: //doi.org/10.1038/s41524-021-00662-x [21] Snoek, J., Ripped, O., Swersky, K., Kiros, R., Satish, N., Sundaram, N., Patwary, M.M.A., Prabhat, Adams, R.P.: Scalable Bayesian optimization using deep neural networks. 32nd International Conference on Machine Learning, ICML 2015 3, 2161–2170 (2015) arXiv:1502.05700 [22] Andrae, A., Edler, T.: On Global Electricity Usage of Communication Technology: Trends to 2030. Challenges 6(1), 117–157 (2015) https://doi.org/10.3390/ challe6010117 [23] Yang, D., Wang, X., Shen, R., Li, Y., Gu, L., Zheng, R., Zhao, J., Tian, X.: Global optimization strategy of prosumer data center system operation based on multiagent deep reinforcement learning. Journal of Building Engineering 91(May), 109519 (2024) https://doi.org/10.1016/j.jobe.2024.109519 [24] Patterson, D., Gonzalez, J., Holzle, U., Le, Q., Liang, C., Munguia, L.M., Rothchild, D., So, D.R., Texier, M., Dean, J.: The Carbon Footprint of Machine Learning Training Will Plateau, Then Shrink. Computer 55(7), 18–28 (2022) https://doi. org/10.1109/MC.2022.3148714 arXiv:2204.05149 [25] Tripp, C.E., Perr-Sauer, J., Gafur, J., Nag, A., Purkayastha, A., Zisman, S., Bensen, E.A.: Measuring the Energy Consumption and Efficiency of Deep Neural Networks: An Empirical Analysis and Design Recommendations. arXiv, 1–25 (2024) arXiv:2403.08151 [26] Vries, A.: The growing energy footprint of artificial intelligence. Joule 7(10), 2191–2194 (2023) https://doi.org/10.1016/j.joule.2023.09.004 [27] Mehta, C., Cherney, M.A., Nellis, S.: Nvidia adds jet fuel to AI optimism with

34

record results, $25 billion buyback. Reuters (2023) [28] Steven, R., Klymenko, O.V., Short, M.: Machine learning-accelerated distributed optimisation methods for optimal power flow: A review. Renewable and Sustainable Energy Reviews 226, 116190 (2026) https://doi.org/10.1016/J.RSER.2025.116190 [29] Mizuyama, H., Morinaga, E., Nonaka, T., Kaihara, T., Cieminski, G., Romero, D. (eds.): Advances in Production Management Systems. Cyber-Physical-Human Production Systems: Human-AI Collaboration and Beyond. IFIP Advances in Information and Communication Technology, vol. 767. Springer, Cham (2026). https://doi.org/10.1007/978-3-032-03542-4 . Conference proceedings [30] Rahman, S., Ahsan, A., Pramanik, N.I.: Climate-Regulating Industrial Ecosystems: An AI-Optimised Framework for Green Infrastructure Performance. Sustainability 17(15), 6891 (2025) https://doi.org/10.3390/su17156891 [31] Naturinda, E., Kemigyisha, F., Gidudu, A., Kabenge, I., Omia, E., Aboth, J.: Quantification of greenhouse gas emissions from livestock using remote sensing & artificial intelligence. Artificial Intelligence in Geosciences 6(2), 100147 (2025) https://doi.org/10.1016/J.AIIG.2025.100147 [32] Xin, Y., Zhu, M., Zhang, H., Wang, X.: High-Entropy Materials: A New Paradigm in the Design of Advanced Batteries vol. 18. Springer, ??? (2026). https://doi. org/10.1007/s40820-025-01842-w [33] Bai, X., Zhang, X.: Artificial Intelligence-Powered Materials Science. Nano-Micro Letters 17(1), 1–30 (2025) https://doi.org/10.1007/s40820-024-01634-8 [34] Liu, Z., Shen, D., Chung, Y.G., Li, W., Li, S.: Accelerated discovery of highperformance MOFs for water adsorption chillers through molecular simulation and machine learning. Chemical Engineering Journal 517, 164419 (2025) https: //doi.org/10.1016/j.cej.2025.164419 [35] Bose, S., Sengupta, D., Rayder, T.M., Wang, X., Kirlikovali, K.O., Sekizkardes, A.K., Islamoglu, T., Farha, O.K.: Challenges and Opportunities: Metal–Organic Frameworks for Direct Air Capture. Advanced Functional Materials 34(43) (2024) https://doi.org/10.1002/adfm.202307478 [36] Luers, A., Koomey, J., Masanet, E., Gaffney, O., Creutzig, F., Lavista Ferres, J., Horvitz, E.: Will AI accelerate or delay the race to net-zero emissions? Nature 628(8009), 718–720 (2024) https://doi.org/10.1038/d41586-024-01137-x [37] Luers, A.: Net zero needs AI — five actions to realize its promise. Nature 644(8078), 871–873 (2025) https://doi.org/10.1038/d41586-025-02641-4 [38] Francom, D., Sansó, B.: BASS: An R package for fitting and performing sensitivity analysis of bayesian adaptive spline surfaces. Journal of Statistical Software 94(8), 1–36 (2020) https://doi.org/10.18637/jss.v094.i08 [39] Duan, T., Avati, A., Ding, D.Y., Thai, K.K., Basu, S., Ng, A., Schuler, A.:

35

NGBoost: Natural gradient boosting for probabilistic prediction. 37th International Conference on Machine Learning, ICML 2020 PartF16814, 2670–2680 (2020) arXiv:1910.03225 [40] Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., Vanderplas, J., Passos, A., Cournapeau, D., Brucher, M., Perrot, M., Duchesnay, E.: Scikit-learn: Machine learning in Python. Journal of Machine Learning Research 12, 2825–2830 (2011) [41] Balandat, M., Karrer, B., Jiang, D.R., Daulton, S., Letham, B., Wilson, A.G., Bakshy, E.: BOTORCH: A framework for efficient Monte-Carlo Bayesian optimization. Advances in Neural Information Processing Systems 2020-December(MC) (2020) arXiv:1910.06403 [42] Hutter, F., Hoos, H.H., Leyton-Brown, K.: Sequential Model-Based Optimization for General Algorithm Configuration. In: Coello, C.A.C. (ed.) Learning and Intelligent Optimization, pp. 507–523. Springer, Berlin, Heidelberg (2011) [43] Styrud, J., Mayr, M., Hellsten, E., Krueger, V., Smith, C.: BeBOP-Combining Reactive Planning and Bayesian Optimization to Solve Robotic Manipulation Tasks. Proceedings - IEEE International Conference on Robotics and Automation, 16459–16466 (2024) https://doi.org/10.1109/ICRA57147.2024.10611468 [44] Denison, D.G.T., Mallick, B.K., Smith, A.F.M.: Bayesian MARS. Statistics and Computing 8(4), 337–346 (1998) https://doi.org/10.1023/A:1008824606259 [45] Li, G., Tian, W., Zhang, H., Chen, B.: Building Energy Models at Different Time Scales Based on Multi-Output Machine Learning. Buildings 12(12) (2022) https://doi.org/10.3390/buildings12122109 [46] Zhang, H., Tian, W., Tan, J., Yin, J., Fu, X.: Sensitivity analysis of multiple timescale building energy using Bayesian adaptive spline surfaces. Applied Energy 363, 123042 (2024) https://doi.org/10.1016/j.apenergy.2024.123042 [47] Ausdemore, M.A., McCombs, A., Ries, D., Zhang, A., Shuler, K., Tucker, J.D., Goode, K., Huerta, J.G.: A probabilistic inverse prediction method for predicting plutonium processing conditions. Frontiers in Nuclear Engineering 1(December), 1–11 (2022) https://doi.org/10.3389/fnuen.2022.1083164 [48] Davis, J.V., Marrs, F.W., Cawkwell, M.J., Manner, V.W.: Machine Learning Models for High Explosive Crystal Density and Performance. Chemistry of Materials 36(22), 11109–11118 (2024) https://doi.org/10.1021/acs.chemmater.4c01978 [49] Chen, H., Liu, J., Shen, G.Q., Feng, Z.: Control of existing tunnel deformation caused by shield adjacent undercrossing construction using interpretable machine learning and multiobjective optimization. Automation in Construction 170, 105943 (2025) https://doi.org/10.1016/j.autcon.2024.105943 [50] Mercado, R., Fu, R.S., Yakutovich, A.V., Talirz, L., Haranczyk, M., Smit, B.: In Silico Design of 2D and 3D Covalent Organic Frameworks for Methane Storage

36

Applications. Chemistry of Materials 30(15), 5069–5086 (2018) https://doi.org/ 10.1021/acs.chemmater.8b01425 [51] Ramakrishnan, R., Dral, P.O., Rupp, M., Von Lilienfeld, O.A.: Quantum chemistry structures and properties of 134 kilo molecules. Scientific Data 1, 1–7 (2014) https://doi.org/10.1038/sdata.2014.22 [52] Luttens, A., Cabeza de Vaca, I., Sparring, L., Brea, J., Martínez, A.L., Kahlous, N.A., Radchenko, D.S., Moroz, Y.S., Loza, M.I., Norinder, U., Carlsson, J.: Rapid traversal of vast chemical space using machine learning-guided docking screens. Nature Computational Science 5(April) (2025) https://doi.org/10.1038/ s43588-025-00777-x [53] Wei, Y., Peng, B., Xie, R., Chen, Y., Qin, Y., Wen, P., Bauer, S., Tung, P.-Y., Raabe, D.: Deep active learning for complex systems. Nature Computational Science 5(September) (2024) https://doi.org/10.1038/s43588-025-00858-x [54] Elmi, M., Pawar, V.M., Shaw, M., Wong, D., Zhan, H., Srinivasan, M.A.: Determining the biomechanics of touch sensation in C. elegans. Scientific Reports 7(1), 12329 (2017) https://doi.org/10.1038/s41598-017-12190-0 [55] Essmann, C.L., Elmi, M., Rekatsinas, C., Chrysochoidis, N., Shaw, M., Pawar, V., Srinivasan, M.A., Vavourakis, V.: The influence of internal pressure and neuromuscular agents on C. elegans biomechanics: an empirical and multi-compartmental in silico modelling study. Frontiers in Bioengineering and Biotechnology 12(March), 1–14 (2024) https://doi.org/10.3389/fbioe.2024.1335788 [56] Rekatsinas, C., Krokidas, P., Vavourakis, V., Essmann, C., Giannakopoulos, G.: A physics-constrained machine learning pipeline for young’s modulus prediction in multi-material hyperelastic cylinders guided by contact mechanics. Advanced Intelligence Discovery (2025). In press [57] Hashin, Z.: Failure Criteria for Unidirectional Fiber Composites. Journal of Applied Mechanics 47(2), 329–334 (1980) https://doi.org/10.1115/1.3153664 [58] Rekatsinas, C., Theodosiou, T., Siorikis, D., Tsiaktanis, K., Chrysochoidis, N., Nastos, C., Saravanos, D.: Micromechanics-based multi-scale framework with strain-rate effects for the simulation of ballistic impact on composite laminates. Journal of Composite Materials 58(27), 2897–2914 (2024) https://doi.org/10. 1177/00219983241283618 [59] Lee, Y.S., Graham, E., Jackson, G., Galindo, A., Adjiman, C.S.: A comparison of the performance of multi-objective optimization methodologies for solvent design. In: Computer Aided Chemical Engineering vol. 46, pp. 37–42. Elsevier, ??? (2019). https://doi.org/10.1016/B978-0-12-818634-3.50007-2 [60] Andrei, E.Y., Efetov, D.K., Jarillo-Herrero, P., et al.: The marvels of moiré materials. Nature Reviews Materials 6, 201–206 (2021) https://doi.org/10.1038/ s41578-021-00284-1

37

[61] Chen, X., Fan, X., Li, L., Zhang, N., Niu, Z., Guo, T., et al.: Moiré engineering of electronic phenomena in correlated oxides. Nature Physics 16(6), 631–635 (2020) [62] Shi, B., Qi, P., Jiang, M., Dai, Y., Lin, F., Zhang, H., Fang, Z.: Exotic physical properties of 2d materials modulated by moiré superlattices. Materials Advances 2(17), 5542–5559 (2021) [63] Papia, E.-M., Kondi, A., Nioras, D., Sofos, F., Constantoudis, V.: Data-driven investigation of analytical relationships in moiré structured materials. Computational Materials Science 261, 114270 (2026) https://doi.org/10.1016/j.commatsci. 2025.114270 [64] Jadaun, P., Soreé, B.: Review of orbital magnetism in graphene-based moire materials. Magnetism 3(3), 245–258 (2023) [65] Pryds, N., Park, D.S., Jespersen, T.S., Yun, S.: Twisted oxide membranes: A perspective. APL Materials 12(1) (2024) [66] Hu, G., Wang, M., Mazor, Y., Qiu, C.W., Alù, A.: Tailoring light with layered and moiré metasurfaces. Trends in Chemistry 3(5), 342–358 (2021) [67] Song, T., Sun, Q.C., Anderson, E., Wang, C., Qian, J., Taniguchi, T., et al.: Direct visualization of magnetic domains and moiré magnetism in twisted 2d magnets. Science 374(6571), 1140–1144 (2021) [68] Chatziathanasiou, G.M., Chrysochoidis, N.A., Rekatsinas, C.S., Saravanos, D.A.: A semi-active shunted piezoelectric tuned-mass-damper for multi-modal vibration control of large flexible structures. Journal of Sound and Vibration 537(July), 117222 (2022) https://doi.org/10.1016/j.jsv.2022.117222 [69] Chatziathanasiou, G.M., Chrysochoidis, N.A., Saravanos, D.A.: A semi-active shunted piezoelectric tuned mass damper for robust vibration control. Journal of Vibration and Control 28(21-22), 2969–2983 (2022) https://doi.org/10.1177/ 10775463211026487 [70] Chatziathanasiou, G.M., Chrysochoidis, N.A., Saravanos, D.A.: Multi-tonal vibration suppression in flexible airframe structures via a semi-active shunted piezoelectric tuned mass damper. Journal of Vibration and Control 0(0) (2025) https://doi.org/10.1177/10775463251341365

38

Frugal Bayesian Optimization: Scalable Surrogates for Data- and Resource-Limited Discovery Panagiotis Krokidas1*, Christoforos Rekatsinas1,2 , Vassilis Sioros1,3 , Grigorios M. Chatziathanasiou4 , Efi-Maria Papia5,6 , George Giannakopoulos1,7 1*

Institute of Informatics and Telecommunications, National Centre for Scientific Research "Demokritos", Agia Paraskevi, Greece. 2 Department of Mechanical Engineering and Aeronautics, University of Patras, Patras Greece. 3 Department of Informatics and Telecommunications, National and Kapodistrian University of Athens, Athens, Greece. 4 School of Mechanical Engineering, Hellenic Mediterranean University, Heraklion, Crete, Greece. 5* Institute of Nanoscience and Nanotechnology, National Centre for Scientific Research "Demokritos", Agia Paraskevi, Greece. 6 Department of Physics, National and Kapodistrian University of Athens, Athens, Greece. 7 SciFY PNPC, Agia Paraskevi, Greece.

*Corresponding author(s). E-mail(s): [email protected]; Contributing authors: [email protected]; [email protected]; [email protected]; [email protected]; [email protected];

1

1 The Benchmarking Functions of this work

Fig. S1 distribution plots for the target property, showing the global optimum and the top-100 values

2

Fig. S2 optimization results in terms of best solution and top-100 solutions identification as a function of sample acquisition

2 AUC areas

Table S1 AUC metrics for the BASS surrogate across benchmark functions.

AUC over time Dataset expschaffer6 rastrigin Michalewicz Ackley Schwefel Styblinski Weierstrass Schaffer7

AUC over samples

Best found

Top–100 acquired

Product

Best found

Top–100 acquired

Product

6.482 31.4 2.486 4.927 2907 58.15 2.013 135

742.2 339.9 313.1 578.5 6.83E+02 210.8 226.8 5.20E+02

4.81E+03 1.07E+04 7.78E+02 2.85E+03 1.99E+06 1.23E+04 4.57E+02 7.02E+04

871.3 4411 366.1 693 4.10E+05 5953 301.2 1.96E+04

9.85E+04 3.82E+04 3.99E+04 5.34E+04 9.28E+04 2.33E+04 2.47E+04 6.63E+04

8.58E+07 1.68E+08 1.46E+07 3.70E+07 3.81E+10 1.38E+08 7.45E+06 1.30E+09

3

Table S2 AUC metrics for the GP surrogate across benchmark functions.

AUC over time Dataset expschaffer6 rastrigin Michalewicz Ackley Schwefel Styblinski Weierstrass Schaffer7

AUC over samples

Best found

Top–100 acquired

Product

Best found

Top–100 acquired

Product

6.278 55.45 3.591 12.18 4333 24 0.7901 24.86

827.4 856.6 788.7 806.1 1.03E+03 592.4 775.4 7.37E+02

5.19E+03 4.75E+04 2.83E+03 9.82E+03 4.47E+06 1.42E+04 6.13E+02 1.83E+04

846.8 9748 559 1650 4.99E+05 5423 348.6 8523

9.17E+04 9.77E+04 3.99E+04 9.14E+04 9.79E+04 7.46E+04 8.69E+04 8.70E+04

7.76E+07 9.52E+08 2.23E+07 1.51E+08 4.89E+10 4.04E+08 3.03E+07 7.42E+08

Table S3 AUC metrics for the NGBoost surrogate across benchmark functions.

AUC over time Dataset expschaffer6 rastrigin Michalewicz Ackley Schwefel Styblinski Weierstrass Schaffer7

AUC over samples

Best found

Top–100 acquired

Product

Best found

Top–100 acquired

Product

0.8508 15.41 0.7943 4.273 1803 63.9 1.7 60.42

188 359.1 188.8 193.7 6.35E+02 341 182.3 529.8

1.60E+02 5.53E+03 1.50E+02 8.28E+02 1.14E+06 2.18E+04 3.10E+02 3.20E+04

96.79 1791 96.63 494.9 2.33E+05 7223 212.7 6788

2.04E+04 3.90E+04 2.14E+04 2.15E+04 8.11E+04 3.77E+04 8.69E+04 5.74E+04

1.97E+06 6.98E+07 2.07E+06 1.06E+07 1.89E+10 2.73E+08 1.85E+07 3.89E+08

Table S4 AUC metrics for the Random Forest surrogate across benchmark functions.

AUC over time Dataset expschaffer6 rastrigin Michalewicz Ackley Schwefel Styblinski Weierstrass Schaffer7

AUC over samples

Best found

Top–100 acquired

Product

Best found

Top–100 acquired

Product

0.6056 16.68 0.6494 2.008 1808 35.63 1.214 36.66

152.1 275.8 194.7 117.9 5.59E+02 195.2 136.3 184.2

9.21E+01 4.60E+03 1.26E+02 2.37E+02 1.01E+06 6.95E+03 1.65E+02 6.75E+03

139.3 3536 158 470.9 3.23E+05 7737 264.1 8403

3.04E+04 5.68E+04 4.03E+04 2.48E+04 9.39E+04 4.18E+04 2.62E+04 4.00E+04

4.23E+06 2.01E+08 6.37E+06 1.17E+07 3.03E+10 3.23E+08 6.92E+06 3.36E+08

4

5

Fig. S3 distribution plots for the target property, showing the global optimum and the top-100 values

Table S5 AUC metrics for the BASS surrogate across real-case datasets.

AUC over time

AUC over samples

Dataset

Best found

Top–100 acquired

Product

Best found

Top–100 acquired

Product

Moire pressure_vessel MOFs proteins lunar worm neural QM9 multi_PZT

0.2773 0.1666 13.15 61.2 187.4 5.353E-04 0.08093 0.05715 44.81

678.5 536.2 289.7 979.2 2372 1010 204.6 297.6 662.9

1.88E+02 8.93E+01 3.81E+03 5.99E+04 4.45E+05 5.41E-01 1.66E+01 1.54E+01 2.97E+04

34.38 29.93 2424 6047 8768 0.05232 16.83 8.26 5591

8.40E+04 9.47E+04 4.98E+04 9.47E+04 9.88E+04 9.90E+04 4.23E+04 4.87E+04 8.16E+04

2.89E+06 2.84E+06 1.21E+08 5.73E+08 8.66E+08 5.18E+03 7.12E+05 4.02E+05 4.56E+08

Table S6 AUC metrics for the GP surrogate across real-case datasets.

AUC over time

AUC over samples

Dataset

Best found

Top–100 acquired

Product

Best found

Top–100 acquired

Product

Moire pressure_vessel MOFs proteins lunar worm neural GM9 multi_PZT

0.08042 0.05793 43.95 754.5 117.9 6.54E-06 0.03758 0.209 43.14

1291 259 1612 1.39E+04 2887 2078 72.85 1719 423.6

1.04E+02 1.50E+01 7.08E+04 1.05E+07 3.40E+05 1.36E-02 2.74E+00 3.59E+02 1.83E+04

38.27 25.46 3024 5996 5360 5.353E-04 12.81 24.73 5828

8.42E+04 9.24E+04 5.41E+04 9.27E+04 9.81E+04 9.26E+04 2.01E+04 7.63E+04 6.32E+04

3.22E+06 2.35E+06 1.64E+08 5.56E+08 5.26E+08 4.96E+01 2.57E+05 1.89E+06 3.69E+08

6

Table S7 AUC metrics for the NGBoost surrogate across real-case datasets.

AUC over time

AUC over samples

Dataset

Best found

Top–100 acquired

Product

Best found

Top–100 acquired

Product

Moire pressure_vessel MOFs proteins lunar worm neural QM9 multi_PZT

0.6535 0.04698 32.03 632.5 181.1 7.51E-05 0.02141 1.505 49.65

1082 219.5 221.8 1.03E+04 4208 1531 29.92 799.2 437.4

7.07E+02 1.03E+01 7.10E+03 6.53E+06 7.62E+05 1.15E-01 6.41E-01 1.20E+03 2.17E+04

52.34 18.52 5268 5787 4197 7814 12.85 81.13 5764

8.79E+04 8.60E+04 3.66E+04 9.39E+04 9.75E+04 9.79E+04 1.81E+04 4.28E+04 5.10E+04

4.60E+06 1.59E+06 1.93E+08 5.43E+08 4.09E+08 7.65E+08 2.32E+05 3.47E+06 2.94E+08

Table S8 AUC metrics for the Random Forest surrogate across real-case datasets.

AUC over time

AUC over samples

Dataset

Best found

Top–100 acquired

Product

Best found

Top–100 acquired

Product

Moire pressure_vessel MOFs proteins lunar worm neural QM9 multi_PZT

0.307 0.03322 4.482 303 99.66 3.31E-05 0.02688 0.08900 19.22

742.2 129.9 99.23 5068 1556 1072 23.37 435.40 418.40

2.28E+02 4.32E+00 4.45E+02 1.54E+06 1.55E+05 3.55E-02 6.28E-01 3.88E+01 8.04E+03

41.63 24.73 2111 5768 7495 5029 43.77 17.20 4336

8.76E+04 9.08E+04 4.38E+04 9.46E+04 9.83E+04 9.85E+04 3.78E+04 6.81E+04 8.96E+04

3.65E+06 2.25E+06 9.24E+07 5.45E+08 7.37E+08 4.96E+08 1.66E+06 1.17E+06 3.89E+08

3 Computational Setup All computations were performed on a personal desktop equipped with an Intel Core i9-10900K CPU, 64 GB of RAM, and an NVIDIA GeForce RTX 3070 Ti GPU (8 GB VRAM), running Windows 11 Pro (64-bit).

7