Paper deep dive
Optimization of time-consuming experimental conditions using pseudo-experimental data guided by adaptive polynomial regression
Hirotaka Sugawara, Yujin Taguchi, Kei Minagawa, Yusuke Hiki, Takashi Morikura, Akira Funahashi
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 90%
Last extracted: 8/2/2026, 1:18:10 PM
Summary
The paper introduces PolyBO, a Bayesian optimization method that accelerates optimization in time-consuming experimental settings by generating high-quality pseudo-experimental data using an adaptively updated polynomial regression model. PolyBO combines real experimental data with these pseudo-data to update the surrogate model, achieving significant reductions in optimization time compared to conventional methods like vanilla BO, TSBO, and BOPP, particularly in high-dimensional spaces and real-world material composition optimization.
Entities (10)
Relation Signals (9)
PolyBO → generates → pseudo-experimental data
confidence 95% · PolyBO performs BO efficiently by generating pseudo-experimental data with an adaptively updated versatile parametric model.
PolyBO → improves → optimization time
confidence 95% · PolyBO reduced the optimization time by a median of 42%... For a real-world material composition optimization problem, PolyBO reduced the optimization time by a median of 96%
PolyBO → uses → Polynomial Regression
confidence 95% · PolyBO performs BO efficiently by generating pseudo-experimental data with an adaptively updated versatile parametric model. This low-capacity polynomial regression model...
PolyBO → comparedwith → vanilla BO
confidence 90% · we found that PolyBO reduced the optimization time by a median of 42%... compared with conventional methods.
PolyBO → updates → Surrogate Model
confidence 90% · PolyBO updates the BO surrogate model with a combined dataset consisting of experimental data and pseudo-experimental data
PolyBO → comparedwith → BOPP
confidence 85% · We compared PolyBO with vanilla BO, random sampling, and BOPP [33]
PolyBO → comparedwith → TSBO
confidence 85% · Algorithms were compared using simple regret... compared with TSBO [47], the latter only in the highest-dimensional setting
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Bayesian optimization (BO) is an optimization method that sequentially proposes the next candidate explainable variables for optimizing target variables by balancing exploration and exploitation. BO is often used under a limited evaluation budget, such as hyperparameter tuning of deep learning. Despite its effectiveness, conventional BO may have poor convergence in practical experimental science where each evaluation is often costly and time-consuming. Recently, BO methods have been proposed that accelerate optimization by using pseudo-experimental data that simulate experimental data. However, when only a limited number of experimental data are available, the generated pseudo-experimental data may be of insufficient quality. In this study, we developed PolyBO to improve optimization time by generating high-quality pseudo-experimental data even when the number of trials is limited. PolyBO performs BO efficiently by generating pseudo-experimental data with an adaptively updated versatile parametric model. This low-capacity polynomial regression model is intended to enable efficient BO even with limited experimental data. PolyBO updates the BO surrogate model with a combined dataset consisting of experimental data and pseudo-experimental data and then performs optimization. Using synthetic benchmark functions with diverse landscapes, we found that PolyBO reduced the optimization time by a median of 42\%. For a real-world material composition optimization problem, PolyBO reduced the optimization time by a median of 96\% compared with conventional methods. Overall, PolyBO achieves efficient optimization in settings where each experiment requires a long time.
Tags
Links
- Source: https://arxiv.org/abs/2607.22238v1
- Canonical: https://arxiv.org/abs/2607.22238v1
Trouble viewing inline? Open PDF directly →
Full Text
80,087 characters extracted from source content.
Expand or collapse full text
Optimization of time-consuming experimental conditions using pseudo-experimental data guided by adaptive polynomial regression Hirotaka Sugawara1 , Yujin Taguchi1 , Kei Minagawa1, Yusuke Hiki1,2 , Takashi Morikura1 , Akira Funahashi1,2,∗ 1 Graduate School of Science and Technology, Keio University, Kanagawa, Japan 2 Department of Biosciences and Informatics, Keio University, Kanagawa, Japan ∗ funa@bio.keio.ac.jp Abstract Bayesian optimization (BO) is an optimization method that sequentially proposes the next candidate explainable variables for optimizing target variables by balancing exploration and exploitation. BO is often used under a limited evaluation budget, such as hyperparameter tuning of deep learning. Despite its effectiveness, conventional BO may have poor convergence in practical experimental science where each evaluation is often costly and time-consuming. Recently, BO methods have been proposed that accelerate optimization by using pseudo-experimental data that simulate experimental data. However, when only a limited number of experimental data are available, the generated pseudo-experimental data may be of insufficient quality. In this study, we developed PolyBO to improve optimization time by generating high-quality pseudo-experimental data even when the number of trials is limited. PolyBO performs BO efficiently by generating pseudo-experimental data with an adaptively updated versatile parametric model. This low-capacity polynomial regression model is intended to enable efficient BO even with limited experimental data. PolyBO updates the BO surrogate model with a combined dataset consisting of experimental data and pseudo-experimental data and then performs optimization. Using synthetic benchmark functions with diverse landscapes, we found that PolyBO reduced the optimization time by a median of 42%. For a real-world material composition optimization problem, PolyBO reduced the optimization time by a median of 96% compared with conventional methods. Overall, PolyBO achieves efficient optimization in settings where each experiment requires a long time. Introduction Optimization of experimental conditions, including material composition, has contributed to the discovery and design of new materials [42, 23, 39]. In particular, computational optimization of the development of biomaterials, which usually requires time-consuming experiments, has accelerated biological insights and applications across various biological fields, ranging from protein engineering [1] to cellular agriculture [34] and biopharmaceutical manufacturing [41, 14]. In the optimization of experimental conditions, the goal is to search for those that maximize the value of the objective function within a search space defined by the explainable variables of the target system. Response surface methodology is a traditional design-of-experiments optimization method that approximates the objective function with a polynomial function and searches for optimal experimental conditions by sampling uniformly across the search space [3]. For example, using this approach, Skrivergaard et al. [36] have optimized the concentrations of three serum-free medium components, fibroblast growth factor 2, fetuin, and bovine serum albumin, and formulated a medium that exceeds the cell-proliferation performance of a conventional medium containing 10% fetal bovine serum. However, (i) accurate modeling of target-system behavior is difficult because nonlinear interactions among explainable variables are intertwined [42, 5, 15], and (i) the relationship between explainable variables and objective values is a black-box relationship, and experiments for evaluating it are time- and resource-intensive [31, 50]. To handle complex nonlinear responses and efficiently identify promising experimental conditions with a limited number of experiments, an alternative method is needed. In biology, Bayesian optimization (BO) [21, 8] meets these requirements [46, 16, 9]. In BO, a surrogate model approximates the objective function, and selection of the next experimental condition is based on its predictive mean and uncertainty. Because BO balances exploration of unobserved conditions and exploitation of past experimental results, it efficiently searches for and proposes experimental conditions [44]. However, vanilla BO faces two challenges when applied to time-consuming experiments such as biological experiments. First, often >10>10 explainable variables are involved [34, 14, 5, 15, 30], and BO is particularly inefficient in such high-dimensional spaces [44, 19, 40]. Latent-space Bayesian optimization (LSBO) is a major strategy for addressing this challenge [40, 10, 11]. LSBO performs BO in a d-dimensional latent space, which is obtained by strongly reducing the dimensionality of a D-dimensional search space. A variational autoencoder (VAE) is used mainly for the dimensionality reduction [20]. Although LSBO mitigates the curse of dimensionality, its search performance depends on whether the latent space appropriately reflects the properties of the search space [40, 10, 11]. If the limited amount of data available for constructing the latent space is insufficient for VAE training, the BO performance is reduced. Second, if a long time is required for each biological experiment, the number of iterations needed for optimization cannot be performed within a practical timeframe. For example, a single testicular tissue culture experiment takes up to 5 weeks [17]. BO often requires dozens of experiments to obtain sufficient improvement; in this case, 50 experiments would take as long as 250 weeks. In short-term experiments (minutes to days), the total time is often acceptable: BO for 30–100 iterations has been performed in experimental systems where each iteration took only 6–20 min [25, 28, 29]. Together, high dimensionality and long duration of experiments can hinder the convergence of optimization. One solution is to augment the dataset with pseudo-experimental data. Yin et al. [47] introduced teacher–student-based semi-supervised learning into LSBO and proposed Teacher–Student Bayesian Optimization (TSBO), in which a teacher model assigns objective values as pseudo-labels to explainable variables sampled in the latent space, and a student model is trained with this information. By correcting the teacher model predictions through feedback from the student model, TSBO improves the quality of pseudo-experimental data and performs LSBO while updating the surrogate model using a combined dataset of experimental and pseudo-experimental data. The effectiveness of TSBO depends on how well the latent representation preserves the structure of the search space because TSBO assumes the LSBO framework. In particular, as mentioned above, insufficient VAE training may degrade BO search performance. This limitation motivates the development of a method that generates pseudo-experimental data in the search space rather than in the latent space. Qian et al. [33] have proposed a BO method (BOPP) that generates pseudo-experimental data around acquired experimental data and updates the surrogate model using a combined dataset of experimental data and pseudo-experimental data. Assuming that nearby explainable variables have close objective values, BOPP uses the objective values of experimental data to assign those to the generated nearby pseudo-experimental data. Because BOPP generates pseudo-experimental data only around experimental data, its effectiveness is limited in regions where the objective function changes drastically. Thus, pseudo-experimental data–based methods still need substantial improvement for application to small experimental datasets. In this study, to improve BO performance with small datasets, we propose PolyBO, a BO method that generates pseudo-experimental data from an adaptively updated versatile parametric model retrained whenever new experimental data are acquired. PolyBO generates pseudo-experimental data in the original search space without relying on latent representations, whose training can become unstable if data are limited. PolyBO assigns objective values to explainable variables sampled from the search space using a polynomial regression model, and it updates the surrogate model with a combined dataset of experimental and pseudo-experimental data. PolyBO design was based on the hypothesis that, in small-dataset settings, polynomial regression shortens optimization time, quantified here as the number of optimization iterations. Polynomial regression is more flexible than a linear model but has lower capacity than deep learning models such as the multilayer perceptron used as the teacher model in TSBO; however, its capacity is sufficient to regress a broad range of functions. To evaluate whether the PolyBO can reduce the optimization time, we used it to optimize 24 synthetic benchmark functions from Black-Box Optimization Benchmarking (BBOB) [13] and compared the time required for PolyBO and vanilla BO to reach the same performance. In high-dimensional setting (D=20D=20), PolyBO demonstrated its effectiveness in a high-dimensional setting by achieving a median 42% reduction in optimization time compared with vanilla BO. To examine whether PolyBO can be applied to a real-world optimization problem, we used it to optimize composition of a material. PolyBO demonstrated its effectiveness in this setting by reaching the performance level achieved by conventional methods in a median 96% shorter optimization time. Materials and Methods Problem setting We consider an optimization problem in which experimental conditions ∈x are varied within the experimentally allowable search space ⊂ℝDX ^D to search for the parameter ∗x^* that maximizes the experimental result, represented as the objective function f:→ℝf:X . ∗=argmaxf()x^*= xargmax\,f(x) (1) PolyBO algorithm PolyBO aims to improve the BO efficiency when each experiment requires a long time; pseudo-experimental data generated from polynomial regression accelerate convergence. In comparison with vanilla BO (Fig. 1a), three steps are added in PolyBO (Fig. 1b). In vanilla BO, a Gaussian process regression model is used as a surrogate model and is trained using experimental data. Using the trained model, an acquisition function is calculated to select the next experimental condition that maximizes the acquisition function, and the next experiment is conducted under the selected condition. The acquisition function proposes the next experimental conditions by balancing exploration of unobserved regions and exploitation of previously observed objective values. Expected Improvement (EI) (Eq. 2) [27] and GP-UCB (Eq. 3) [37] have been proposed as representative acquisition functions. EI calculates the expected improvement over the known best value f(+)f(x^+), whereas GP-UCB calculates an upper confidence bound from the posterior mean μ()μ(x) and standard deviation σ()σ(x) of the trained Gaussian process regression. αEI()=(μ()−f(+))Φ(Z)+σ()ϕ(Z)σ()>00σ()=0,Z=μ()−f(+)σ() _EI(x)= cases(μ(x)-f(x^+)) (Z)+σ(x)φ(Z)&σ(x)>0\\ 0&σ(x)=0 cases, Z= μ(x)-f(x^+)σ(x) (2) αGP-UCB()=μ()+βσ() _GP-UCB(x)=μ(x)+ βσ(x) (3) Here, α denotes the acquisition function; ϕφ and Φ are the probability density function and cumulative distribution function, respectively, of the standard normal distribution; Z is the standardized variable, and β is a hyperparameter that controls the exploration–exploitation balance. In PolyBO, a polynomial regression model is trained using experimental data, pseudo-experimental data are generated by sampling from the trained polynomial regression model, and the Gaussian process regression model is trained using the combined experimental and pseudo-experimental data, which are treated equally in the same format as pairs of explainable variables and objective values. The acquisition function is then calculated as in vanilla BO, and the next experiment is conducted under the condition that maximizes it. The pseudocode of the PolyBO algorithm is shown in Algorithm 1. At iteration t=0t=0, an initial experimental dataset of size k, 0:0=(i(t=0),yi(t=0))i=1kD_0:0=\(x^(t=0)_i,y^(t=0)_i)\_i=1^k, is prepared for training the Gaussian process regression model. Here, k is the initial dataset size, and the index i=1,…,ki=1,…,k enumerates the samples in the initial dataset. The experimental data consist of explainable variables i(0)x^(0)_i randomly sampled from the search space and the corresponding objective values yi(0)y^(0)_i. The set consisting of the initial data and the experimental data obtained up to iteration t−1t-1 is denoted as 0:t−1=0:0∪((τ),y(τ))τ=1t−1D_0:t-1=D_0:0∪\(x^(τ),y^(τ))\_τ=1^t-1. Here, τ=1,…,t−1τ=1,…,t-1 indexes the experimental data sequentially acquired after the initial dataset. When t=1t=1, ((τ),y(τ))τ=1t−1\(x^(τ),y^(τ))\_τ=1^t-1 is the empty set, and therefore 0:t−1=0:0D_0:t-1=D_0:0. The pthp^th-degree polynomial regression model fp()f_p(x) is trained using the experimental data 0:t−1D_0:t-1. Polynomial regression expands the explainable variables x into a new feature vector ϕ()φ(x) and then applies linear regression; ϕ()φ(x) includes not only power terms up to degree p for each explainable variable but also interaction terms, expressed as products of different explainable variables. Let be the design matrix constructed by stacking ϕ()⊤φ(x) for each training sample and =[y1(0),…,yk(0),y(1),…,y(t−1)]⊤y=[y_1^(0),…,y_k^(0),y^(1),…,y^(t-1)] be the corresponding objective values; the model weight parameter w is then analytically obtained by solving the normal equation (Eq. 4). =(⊤)−1⊤w=( )^-1 y (4) The number of parameters in a model that includes all power and interaction terms is (D+p)=(D+1)(D+2)⋯(D+p)p!. D+pp= (D+1)(D+2)·s(D+p)p!. Pseudo-experimental data t′=(′j(t),yj′(t))j=1m′D_t =\(x ^(t)_j,y (t)_j)\_j=1^m are generated using the trained model fp()f_p(x). Here, m′m is the pseudo-experimental dataset size, and j is the sample index. The pseudo-experimental data consist of explainable variables ′j(t)x ^(t)_j randomly sampled from the search space and objective values yj′(t)=fp(′j(t))y (t)_j=f_p(x ^(t)_j) calculated by substituting them into the trained model fpf_p. The posterior distribution is predicted by Gaussian process regression using 0:t−1∪t′D_0:t-1 _t , which combines the pseudo-experimental data t′D_t and experimental data 0:t−1D_0:t-1, and then the acquisition function is calculated. When using EI as the acquisition function, PolyBO uses pseudo-experimental data only as auxiliary information for constructing the Gaussian process posterior and defines the incumbent value f(+)f(x^+) in Eq. 2 solely from the experimental data 0:t−1D_0:t-1. This design prevents pseudo-experimental objective values potentially affected by regression error from directly determining the EI improvement baseline. A limited-memory Broyden–Fletcher–Goldfarb–Shanno algorithm with bound constraints [4] was used to numerically find the maximum of the acquisition function. The explainable variable (t)x^(t) that maximizes the acquisition function is obtained, and then the experiment is conducted. Using the experimental result y(t)y^(t) as the objective value, the experimental data are updated as 0:t=0:t−1∪((t),y(t))D_0:t=D_0:t-1∪\(x^(t),y^(t))\, and then the iteration proceeds. PolyBO optimizes experimental conditions by repeating the above process. At every iteration, the update mechanism in PolyBO replaces the pseudo-experimental dataset with a newly generated dataset t′D_t of size m′m and discards pseudo-experimental data generated in the previous iterations. This prevents low-quality pseudo-experimental data generated when only a few experimental data are available from degrading convergence performance in the middle and later stages of optimization. Algorithm 1 PolyBO algorithm 1: Set initial dataset 0:0=(i(0),yi(0))i=1kD_0:0=\(x^(0)_i,y^(0)_i)\_i=1^k, where i(0)x^(0)_i are randomly sampled from the search space and yi(0)y^(0)_i are obtained from experiments. 2: for t=1,2,…t=1,2,… do 3: Train pthp^th-degree polynomial regression model fp()f_p(x) with dataset 0:t−1D_0:t-1. 4: Sample m′m explainable variables ′j(t)j=1m′\x ^(t)_j\_j=1^m randomly from the search space. 5: for j=1,2,…,m′j=1,2,…,m do 6: Compute yj′(t)=fp(′j(t))y (t)_j=f_p(x ^(t)_j). 7: end for 8: Set t′=(′j(t),yj′(t))j=1m′D_t =\(x ^(t)_j,y (t)_j)\_j=1^m . 9: Update the Gaussian process regression model with 0:t−1∪t′D_0:t-1 _t . 10: Compute the acquisition function α()α(x) using the updated Gaussian process regression model. 11: Select the next point (t)=argmaxα()x^(t)=argmax_xα(x). 12: Conduct an experiment at (t)x^(t) to obtain y(t)y^(t). 13: Update the dataset 0:t=0:t−1∪((t),y(t))D_0:t=D_0:t-1∪\(x^(t),y^(t))\. 14: end for Numerical experiments using synthetic benchmark functions To evaluate the performance of PolyBO, we optimized 24 BBOB functions [13]. These functions are continuous, with diverse landscapes defined on [−5,5]D[-5,5]^D for arbitrary dimensionality D, and they can be regarded as broadly simulating biological phenomena with high-dimensional and complex responses. These functions were implemented using COmparing Continuous Optimizers (COCO) [12]; we selected D=2D=2 and 55 as low-dimensional settings and D=10D=10 and 2020 as high-dimensional settings. For each function, we searched for optimal conditions using different initial conditions generated from random numbers based on the Philox algorithm [35] with 30 random seeds and evaluated the optimization process. Because PolyBO aims to improve optimization efficiency when the amount of data is small, the initial dataset size was set to two. Algorithms were compared using simple regret (Eq. 5), defined as the difference between the instance-specific reference optimal value f(∗)f(x^*), predefined by the benchmark, and the maximum value found after n searches. simple regret=f(∗)−maxi=1,⋯,nf(i)simple regret=f(x^*)- i=1,\;·s,\;nmaxf(x^i) (5) A smaller simple regret indicates a more optimal value. EI (Eq. 2) was used as the acquisition function. All methods were implemented on the basis of BoTorch [2], and the source code and data are available at https://github.com/funalab/PolyBO. Performance comparison with vanilla BO We compared PolyBO with vanilla BO, random sampling, and BOPP [33], as well as with TSBO [47], the latter only in the highest-dimensional setting (D=20D=20). The maximum number of all searches from the initial conditions was set to 100. We used the evaluation metric I50I_50 [38], which is the number of iterations required for the compared method to reach the simple regret reached by vanilla BO after 50 optimization iterations. I50<50I_50<50 indicates a shorter optimization time than that of vanilla BO, and small I50I_50 indicates high optimization performance. Analysis of the optimization process To compare the optimization process between vanilla BO and PolyBO, the sampled points at each iteration were visualized in a reduced-dimensional space. As dimensionality-reduction algorithms, we used t-distributed Stochastic Neighbor Embedding (t-SNE) [24], which emphasizes the preservation of the local neighborhood structure, and Uniform Manifold Approximation and Projection (UMAP) [26], which tends to preserve global structures such as inter-cluster distances more accurately than t-SNE. To sufficiently cover the entire high-dimensional search space, 40,000 explainable variables were obtained by Latin hypercube sampling [22], and dimensionality was reduced on these variables together with all sampled points from vanilla BO and PolyBO. We also quantified optimization behavior to support the dimensionality reduction–based assessment. As evaluation metrics, we used Observation Entropy (OE) (Eq. 6) [32] and Proposed-point Entropy (PE) (Eq. 7), which we devised, and compared them between vanilla BO and PolyBO. OEt=Dt∑i=1tlogϵik+ψ(t)−ψ(1)+logVDOE_t= Dt _i=1^t _i^k+ψ(t)-ψ(1)+ V_D (6) PEt=−∫p(y∣(t),0:t−1∪t′)logp(y∣(t),0:t−1∪t′)y=12log(2πeσt2((t))).PE_t=- p(y ^(t),D_0:t-1 _t ) p(y ^(t),D_0:t-1 _t )\,dy= 12 (2π e\, _t^2(x^(t)) ). (7) Here, ϵik _i^k is the distance from point (i)x^(i) to its kthk^th nearest neighbor, VD=πD/2/Γ(1+D/2)V_D=π^D/2/ (1+D/2) is the volume of a D-dimensional unit ball, Γ is the gamma function, and ψ is the digamma function. Following Papenmeier et al. [32], we set k=log(t)k= (t). OE quantifies the density of explainable variables among the sampled points obtained up to iteration t, whereas PE quantifies the predictive variance at the explainable variable proposed at iteration t. In BO, it is important to consider both explainable variables and objective values. OE focuses on explainable variables and PE on objective values, i.e. these metrics quantify complementary aspects of uncertainty. For both metrics, higher values indicate a more exploration-oriented method. Robustness analysis of hyperparameters and acquisition functions The PolyBO-specific hyperparameters are the pseudo-experimental dataset of size m′m and the maximum polynomial degree p. We assessed how changes in the hyperparameters affect PolyBO performance (sensitivity analysis). At D=20D=20, one hyperparameter was fixed, while the other was varied. We measured I50I_50 at p=4p=4 and m′m varied over 1, 5, 10, 15, 20, 25, 40, 80, 100, 200\1,\;5,\;10,\;15,\;20,\;25,\;40,\;80,\;100,\;200\, and at m′=10m =10 and p varied over 1, 2, 3, 4, 5, 6, 7, 8\1,\;2,\;3,\;4,\;5,\;6,\;7,\;8\. To examine whether the optimization performance of PolyBO depends on the EI acquisition function (Eq. 2), we compared the optimization process with that obtained using GP-UCB (Eq. 3). The hyperparameter β affects the optimization strategy and was fixed at β=1β=1 as in [7, 49]; it was also used by Papenmeier et al. [32] because β has no marked bias toward exploratory or exploitative behavior. Ablation study of the update mechanism The update mechanism of PolyBO (i) discards past pseudo-experimental data and (i) maintains a constant pseudo-experimental dataset size. To clarify the effect of the former on optimization performance, we evaluated how a mechanism that does not discard pseudo-experimental data and uses all such previously generated data (the without-reset mechanism) affects optimization performance. To clarify the effect of the latter, we evaluated the m′m -scaling mechanism, which does not retain pseudo-experimental data from previous iterations but instead generates a fresh pseudo-experimental dataset of size t×m′t× m at iteration t. Application to a real-world material composition optimization problem Simulated material composition results generated by a neural-network predictor for high-entropy alloys (alloys formed by mixing five or more elements in similar proportions; HEAs) [45, 43] were used as a proxy for experimental results. This simulation model maps composition parameters of up to 10 elements to key mechanical properties: yield strength (σY _Y), ultimate tensile strength (σU _U), and elongation (ε ). Although HEAs are not the biomaterials that are the main target of this study, we selected them as a representative example of an expensive-to-evaluate optimization problem because, like biomaterials, their synthesis requires a long time. We optimized composition parameters for 10 elements. For evaluation, we used the HEA figure of merit (FOMHEA FOM_HEA) (Eq. 8). FOMHEA FOM_HEA represents the inherent strength–ductility trade-off as a weighted combination of the three normalized properties. As in Xian et al. [45], we evaluated FOMHEA FOM_HEA against Experimental Iteration, which is the number of actual experiments performed. FOMHEA=13(σYσYN+σUσUN+εN) FOM_ HEA= 13 ( σ_Yσ_Y_N+ σ_Uσ_U_N+ _N ) (8) The subscript N indicates the normalized value of the corresponding property. We compared PolyBO with BO-EI and the reinforcement-learning method RL-DQN [45]. A search for optimal conditions using different initial conditions was based on 30 random seeds (the HEA N predictor for PolyBO and BO-EI) or 24 random seeds (RL-DQN). We used the same 20 initial conditions for PolyBO as for BO-EI. The upper limit of the number of searches was set to 1500, including the initial data. EI (Eq. 2) was used as the acquisition function. Because this was a discrete and constrained optimization problem, we extended PolyBO to handle this search space (Supplementary Algorithm 1). Following BO-EI, we developed an algorithm in which the atomic fraction xix_i of each element is restricted to a discrete grid with step size δ and simultaneously satisfies the equality constraint ∑i=1dxi=1 _i=1^dx_i=1 as well as the lower and upper bound constraints li≤xi≤uil_i≤ x_i≤ u_i for each element. We generated pseudo-experimental data under the constraint that the sampled explainable variables were drawn uniformly from the discrete grid δG_δ, which satisfied the above constraints. This sampling scheme generated pseudo-experimental data in the same space as experimental data and was expected to provide a more appropriate training-data distribution for Gaussian process regression than sampling from the corresponding continuous constrained space. We quantified PolyBO performance using the number of iterations it required to reach the FOMHEA FOM_HEA attained by RL-DQN after 1500 optimization iterations (I1500I_1500). Results Performance of PolyBO on synthetic benchmark functions To evaluate PolyBO, we compared it with conventional BO, random sampling, and BOPP, which generates pseudo-experimental data around acquired experimental data [33]. At D=20D=20, we also compared PolyBO with TSBO, which incorporates semi-supervised learning [47]. Across the 24 benchmark functions, PolyBO tended to achieve lower simple regret earlier than vanilla BO at D=10D=10 and D=20D=20, but not at D=2D=2 or D=5D=5 (Fig. 2, Supplementary Figs. 1–4). At D=20D=20, PolyBO decreased I50I_50 for all benchmark functions except F05, F13, F16, F21, and F23. This indicated that PolyBO required fewer iterations for optimization (Table 1). The median rate of increase was 12% at D=2D=2, whereas the median rates of reduction were 6% at D=5D=5, 31% at D=10D=10, and 42% at D=20D=20. At D=5D=5, the value could not be calculated for any benchmark functions other than F04, F05, F06, and F16 because the reference simple regret was not reached within 100 iterations. These results indicated that the optimization performance of PolyBO was low in low-dimensional settings, particularly at D=5D=5, but was high in high-dimensional settings. Analysis of the optimization process We visualized optimization processes for different functions at various dimensionalities using scatter plots. As representative examples, the optimization processes for F04 and F06 (D=20D=20) are shown in Figure 3 (t-SNE) and Supplementary Figures 5 (UMAP), 6 (D=2D=2), and 7 (D=5D=5). The entire search spaces for F04 and F06 are shown in Supplementary Fig. 8. At D=20D=20, the tendency of PolyBO to explore a broader region of the search space than vanilla BO (Fig. 3 and Supplementary Fig. 5) indicated that PolyBO is an exploration-oriented method. The t-SNE and UMAP results showed little difference in the global trends, supporting this interpretation. At D=2D=2 and D=5D=5, little difference was found in the extent of the search space covered by PolyBO and vanilla BO (Supplementary Figs. 6 and 7). At D=20D=20, PolyBO found promising regions near the global optimum between iterations 20 and 40 for F04 and between iterations 0 and 20 for F06. We quantified optimization behavior using OE for uncertainty in the input space and PE for predictive uncertainty in objective values (Fig. 4, Supplementary Figs. 9 and 10). At D=2D=2, 55, and 2020, both metrics showed that PolyBO was more exploration-oriented than vanilla BO. These quantitative results supported the qualitative results from the search-space visualization for high-dimensional settings but not for low-dimensional settings. Robustness analysis of hyperparameters and acquisition functions To evaluate the robustness of the PolyBO-specific hyperparameters m′m and p, we assessed the effects of their variation on PolyBO performance through sensitivity analysis (Tables 2 and 3). To examine whether the optimization performance of PolyBO depends on the EI acquisition function, we compared it with that obtained using GP-UCB (Fig. 5). At a maximum of 100 iterations, PolyBO had better optimization performance than vanilla BO over the tested range of pseudo-experimental dataset sizes up to 200 (Table 2). This suggests that PolyBO is unlikely to underperform vanilla BO even when the pseudo-experimental dataset size is increased. Optimization was maximized at m′=5–25m =5--25. This suggests a nonlinear relationship between optimization performance and m′m . At p≠1p≠ 1, PolyBO had better optimization performance than vanilla BO (Table 3). The median optimization performance, ranked from best to worst, was achieved at p=4, 8, 6, 7, 2, 5, 3, 1p=4,\;8,\;6,\;7,\;2,\;5,\;3,\;1. The suggestion was that even p values (2, 4, 6, 8)(2,\;4,\;6,\;8) tended to result in higher optimization performance than the odd values (1, 3, 5, 7)(1,\;3,\;5,\;7). For benchmark functions other than F05, F21, and F22, vanilla BO (EI) and vanilla BO (GP-UCB) showed similar trends, as did PolyBO (EI) and PolyBO (GP-UCB) (Fig. 5). For F05, F21, and F22, vanilla BO (GP-UCB) outperformed vanilla BO (EI), and PolyBO (GP-UCB) outperformed PolyBO (EI); PolyBO (GP-UCB) optimization performance was comparable to or higher than that of vanilla BO (GP-UCB). These results indicate that the performance improvement of PolyBO depends little on differences in acquisition functions. Ablation study of the update mechanism We analyzed the two components of the update mechanism, namely discarding pseudo-experimental data and maintaining constant m′m (Fig. 6). Optimization performance was higher in the model with the update mechanism than in those with the without-reset mechanism or the m′m -scaling mechanism, which does not maintain constant m′m ; the difference was greater in the former than in the latter case. These results show that both discarding pseudo-experimental data and maintaining constant m′m are important and contribute to performance improvement, especially through the update mechanism that discards pseudo-experimental data generated in previous iterations and updates the pseudo-experimental dataset at each iteration. Application to a real-world material composition optimization problem We optimized material composition and compared the optimization processes (Fig. 7). The optimization time was shorter for PolyBO than for the conventional methods. The median value of the metric I1500I_1500 was 59. This indicated that PolyBO reached the performance level achieved by RL-DQN with a 96% shorter optimization time. This result suggests that PolyBO can accelerate optimization when applied not only to synthetic benchmark functions but also to experimentally grounded, real-world search spaces. Discussion We developed PolyBO, an optimization method that uses regression-based pseudo-experimental data to reduce the optimization time. In PolyBO, BO is conducted with augmented data from an adaptively updated versatile parametric model, and these pseudo-experimental data accelerate convergence even when each experiment requires a long time (Fig. 1). In high-dimensional search spaces, PolyBO found more optimal parameters in fewer iterations than did conventional methods (Figs. 2 and 7) and reached comparable conditions in fewer iterations (Table 1). These results demonstrate that PolyBO improves optimization performance in high-dimensional settings. Although this study did not directly evaluate biomaterials, PolyBO may also be applicable to time-consuming biological experiments. By reducing the optimization time, PolyBO can enable searches for experimental conditions that were previously difficult to conduct in time-consuming experimental systems, such as testicular tissue culture, where a single experiment takes up to 5 weeks [17]. Visualization of the optimization process showed that PolyBO found promising regions near the global optimum by emphasizing exploration in the early phase (Fig. 3). The high performance of PolyBO on the F04 and F06 functions (Table 1, Supplementary Fig. 4) can be attributed to its exploration-oriented behavior that allowed it to find promising regions earlier than did vanilla BO. Thus, in high-dimensional settings with broad search spaces, the stronger exploratory behavior of PolyBO may have improved optimization performance by enabling earlier discovery of promising regions. In low-dimensional settings, the discrepancy between the qualitative results from scatter-plot visualization of the optimization process (Fig. 3) and the quantitative results from OE and PE (Fig. 4) arises from the different aspects of optimization behavior: the scatter-plots show the search range of sampled points, whereas OE evaluates density from inter-point distances, and PE evaluates the predictive uncertainty of proposed points. In low-dimensional settings, the search space is relatively small, so vanilla BO can cover a broad search range with fewer sampled points than in high-dimensional settings. As a result, differences between PolyBO and vanilla BO are not obvious in the scatter-plots. However, even when the search range of sampled points is comparable between PolyBO and vanilla BO, PolyBO may select sparse regions far from existing sampled points or points with high predictive uncertainty more frequently than vanilla BO, and these differences may have been reflected in OE and PE. Sensitivity analysis provided guidance for selecting the maximum polynomial degree p. When p=1p=1, the convergence performance of PolyBO was low (Table 3). This was likely because the polynomial regression model has only 21 parameters and low model expressiveness when p=1p=1, preventing the regression model from capturing complex benchmark-function landscapes. When p was varied from 1 to 8, the only benchmark function for which p=1p=1 showed the highest optimization performance was F05, whose search space is mostly linear (Table 3). A linear regression model (p=1p=1) is advantageous when the response of the objective function to the explainable variables is expected to be linear, whereas nonlinear regression models (p≠1p≠ 1) are advantageous for nonlinear responses. For F08, F09, F10, and F20, optimization performance was equal to or higher at even degrees (p=2, 4, 6, 8p=2,\;4,\;6,\;8) than at odd degrees (p=1, 3, 5, 7p=1,\;3,\;5,\;7) (Table 3). F08 and F09 have bent valleys, F10 has a curvature structure with strong ill-conditioning, F20 has a complex global structure near the boundary [13], and pseudo-experimental data generated by even-degree polynomials may be more advantageous for optimization than asymmetric pseudo-experimental data that can arise from odd-degree polynomials. In the ablation study, optimization performance was lower with the without-reset mechanism and with the m′m -scaling mechanism than with the update mechanism. The data in Table 2 supported this result: optimization performance was maximized when the pseudo-experimental dataset size was 5–25; because the m′m -scaling mechanism increases this size to t×m′t× m at iteration t, it exceeds this favorable range after the early iterations when m′=10m =10, explaining the reduced optimization performance shown in Fig. 6. In low-dimensional settings, particularly at D=5D=5, PolyBO optimization performance was low. The suggestion was that it depends on the balance between the size of the search space and the regression error of the pseudo-experimental data. At D=2D=2, the exploration-enhancing property of PolyBO may not have contributed substantially because the search space was small, and vanilla BO could easily reach promising regions. In low-dimensional settings, the expectation that the regression error of the polynomial regression model will be relatively small makes it less likely that pseudo-experimental data strongly impede optimization; consequently, PolyBO performance remained comparable to that of vanilla BO or was lower only in limited cases. At D=5D=5, vanilla BO may already have sufficient convergence performance, whereas its broader search space than at D=2D=2 makes the influence of pseudo-experimental data non-negligible. To evaluate the memory complexity of the PolyBO algorithm, we measured peak memory usage required to solve the normal equation (Eq. 4) while increasing the maximum degree p from 1 to 9 at D=20D=20, and we observed exponential growth (Supplementary Fig. 11a). As the maximum degree p increases, the number of parameters increases as O(D+p)O D+pp (Supplementary Fig. 11b), and the memory usage for solving the normal equation (Eq. 4) in polynomial regression theoretically increases as O((D+p)2)O ( D+pp^2 ). Although the optimal maximum degree differs among benchmark functions (Table 3), increasing the maximum degree p exponentially increases memory complexity, and the optimal degree must therefore be determined while considering hardware specifications. Our data indicate that an optimal pseudo-experimental dataset size rather than addition of a large amount of pseudo-experimental data (Table 2) and its sequential updating are important for maximizing optimization performance (Fig. 6). It is thus important to incorporate the low-fidelity source provided by the polynomial regression model into the surrogate model with an appropriate strength and to update this source as iterations proceed. This approach is similar to that used in multi-fidelity Bayesian optimization (MFBO) [18, 48, 6], which improves the efficiency of optimizing a high-fidelity objective function by using multiple information sources with different evaluation costs and accuracies. If the pseudo-experimental dataset size is too small, enhancement of exploration is insufficient; if it is too large, there is a risk that the surrogate model becomes biased by over-reliance on pseudo-experimental data. In future work, it will be important to extend PolyBO to a framework that sequentially controls the influence of pseudo-experimental data according to regression uncertainty and consistency with experimental data, rather than fixing the size of the pseudo-experimental dataset. Because the fidelity of pseudo-experimental data in PolyBO changes as iterations proceed and many existing MFBO methods assume fixed fidelity for each information source [18, 48, 6], extending MFBO to model iteration-dependent fidelity will be an important direction for future work. References [1] F. H. Arnold (2017) Directed evolution: bringing new chemistry to life. Angewandte Chemie (International Ed. in English) 57 (16), p. 4143. Cited by: Introduction. [2] M. Balandat, B. Karrer, D. R. Jiang, S. Daulton, B. Letham, A. G. Wilson, and E. Bakshy (2020) BoTorch: A Framework for Efficient Monte-Carlo Bayesian Optimization. In Advances in Neural Information Processing Systems 33, External Links: Link Cited by: Numerical experiments using synthetic benchmark functions. [3] M. A. Bezerra, R. E. Santelli, E. P. Oliveira, L. S. Villar, and L. A. Escaleira (2008) Response surface methodology (RSM) as a tool for optimization in analytical chemistry. Talanta 76 (5), p. 965–977. External Links: Document Cited by: Introduction. [4] R. H. Byrd, P. Lu, J. Nocedal, and C. Zhu (1995) A limited memory algorithm for bound constrained optimization. SIAM Journal on scientific computing 16 (5), p. 1190–1208. Cited by: PolyBO algorithm. [5] Z. Cosenza, R. Astudillo, P. I. Frazier, K. Baar, and D. E. Block (2022) Multi-information source Bayesian optimization of culture media for cellular agriculture. Biotechnology and Bioengineering 119 (9), p. 2447–2458. External Links: ISSN 1097-0290, Link, Document Cited by: Introduction, Introduction. [6] B. Do and R. Zhang (2025) Multifidelity Bayesian Optimization: A Review. AIAA Journal 63 (6), p. 2286–2322. External Links: Document, Link, https://doi.org/10.2514/1.J063812 Cited by: Discussion. [7] J. Fromer, R. Wang, M. Manjrekar, A. Tripp, J. M. Hernández-Lobato, and C. W. Coley (2025) Batched Bayesian optimization with correlated candidate uncertainties. External Links: Link Cited by: Robustness analysis of hyperparameters and acquisition functions. [8] R. Garnett (2023) Bayesian optimization. Cambridge University Press. Cited by: Introduction. [9] F. Gisperg, R. Klausser, M. Elshazly, J. Kopp, E. P. Brichtová, and O. Spadiut (2025) Bayesian Optimization in Bioprocess Engineering—Where Do We Stand Today?. Biotechnology and Bioengineering 122 (6), p. 1313–1325. Cited by: Introduction. [10] R. Gómez-Bombarelli, J. N. Wei, D. Duvenaud, J. M. Hernández-Lobato, B. Sánchez-Lengeling, D. Sheberla, J. Aguilera-Iparraguirre, T. D. Hirzel, R. P. Adams, and A. Aspuru-Guzik (2018) Automatic Chemical Design Using a Data-Driven Continuous Representation of Molecules. ACS Central Science 4 (2), p. 268–276. External Links: ISSN 2374-7943, Document, Link Cited by: Introduction. [11] A. Grosnit, R. Tutunov, A. M. Maraval, R. Griffiths, A. I. Cowen-Rivers, L. Yang, L. Zhu, W. Lyu, Z. Chen, J. Wang, et al. (2021) High-dimensional Bayesian optimisation with variational autoencoders and deep metric learning. arXiv preprint arXiv:2106.03609. Cited by: Introduction. [12] N. Hansen, A. Auger, R. Ros, O. Mersmann, T. Tušar, and D. Brockhoff (2021) COCO: A Platform for Comparing Continuous Optimizers in a Black-Box Setting. Optim. Methods and Softw. 36, p. 114–144. External Links: Document Cited by: Numerical experiments using synthetic benchmark functions. [13] N. Hansen, S. Finck, R. Ros, and A. Auger (2009) Real-Parameter Black-Box Optimization Benchmarking 2009: Noiseless Functions Definitions. [Research Report] R-6829, INRIA. inria-00362633v2, 16p.. Cited by: Introduction, Numerical experiments using synthetic benchmark functions, Discussion. [14] T. Hashizume, Y. Ozawa, and B. Ying (2023) Employing active learning in the optimization of culture medium for mammalian cells. npj Systems Biology and Applications 9 (1), p. 20. External Links: ISSN 2056-7189, Document, Link Cited by: Introduction, Introduction. [15] T. Hashizume and B. Ying (2024) Challenges in developing cell culture media using machine learning. Biotechnology Advances 70, p. 108293. External Links: ISSN 0734-9750, Document, Link Cited by: Introduction, Introduction. [16] R. Hu, L. Fu, Y. Chen, J. Chen, Y. Qiao, and T. Si (2022) Protein Engineering via Bayesian Optimization-Guided Evolutionary Algorithm and Robotic Experiments. Briefings in Bioinformatics 24 (1), p. bbac570. External Links: https://academic.oup.com/bib/article-pdf/24/1/bbac570/48782188/bbac570.pdf, ISSN 1477-4054, Document Cited by: Introduction. [17] M. Kamoshita, H. Shirai, H. Nakamura, T. Kishimoto, Y. Hatanaka, D. Mashiko, K. Esashika, J. Yang, S. Yamasaki, T. Ogawa, H. Kimura, and M. Ikawa (2025-01-03) Development of the membrane ceiling method for in vitro spermatogenesis. Scientific Reports 15 (1), p. 625. External Links: ISSN 2045-2322, Document, Link Cited by: Introduction, Discussion. [18] K. Kandasamy, G. Dasarathy, J. Oliva, J. Schneider, and B. Poczos (2019) Multi-fidelity gaussian process bandit optimisation. Journal of Artificial Intelligence Research 66, p. 151–196. Cited by: Discussion. [19] K. Kandasamy, J. Schneider, and B. Poczos (2015) High Dimensional Bayesian Optimisation and Bandits via Additive Models. In Proceedings of the 32nd International Conference on Machine Learning, F. Bach and D. Blei (Eds.), Proceedings of Machine Learning Research, Vol. 37, Lille, France, p. 295–304. External Links: Link Cited by: Introduction. [20] D. P. Kingma and M. Welling (2013) Auto-encoding variational bayes. arXiv preprint arXiv:1312.6114. Cited by: Introduction. [21] H. J. Kushner (1964) A New Method of Locating the Maximum Point of an Arbitrary Multipeak Curve in the Presence of Noise. J. Basic Eng. 86 (1), p. 97–106. External Links: Document Cited by: Introduction. [22] H. Liang, F. Zhao, T. Xu, and J. Zhang (2025) A Deep Learning-Driven Black-Box Benchmark Generation Method via Exploratory Landscape Analysis. Applied Sciences 15 (15). External Links: Link, ISSN 2076-3417, Document Cited by: Analysis of the optimization process. [23] T. Lookman, P. V. Balachandran, D. Xue, and R. Yuan (2019) Active learning in materials science with emphasis on adaptive sampling using uncertainties for targeted design. npj Computational Materials 5 (1), p. 21. External Links: Document, Link Cited by: Introduction. [24] L. v. d. Maaten and G. Hinton (2008) Visualizing data using t-SNE. Journal of machine learning research 9 (Nov), p. 2579–2605. Cited by: Analysis of the optimization process. [25] B. P. MacLeod, F. G. L. Parlane, T. D. Morrissey, F. Häse, L. M. Roch, K. E. Dettelbach, R. Moreira, L. P. E. Yunker, M. B. Rooney, J. R. Deeth, V. Lai, G. J. Ng, H. Situ, R. H. Zhang, M. S. Elliott, T. H. Haley, D. J. Dvorak, A. Aspuru-Guzik, J. E. Hein, and C. P. Berlinguette (2020) Self-driving laboratory for accelerated discovery of thin-film materials. Science Advances 6 (20), p. eaaz8867. External Links: Document, Link, https://w.science.org/doi/pdf/10.1126/sciadv.aaz8867 Cited by: Introduction. [26] L. McInnes, J. Healy, and J. Melville (2018) Umap: Uniform manifold approximation and projection for dimension reduction. arXiv preprint arXiv:1802.03426. Cited by: Analysis of the optimization process. [27] J. Močkus (1975) On Bayesian Methods for Seeking the Extremum. In Optimization Techniques IFIP Technical Conference: Novosibirsk, July 1–7, 1974, p. 400–404. External Links: ISBN 978-3-662-38527-2, Document, Link Cited by: PolyBO algorithm. [28] N. Munyebvu, S. Dunn, and P. D. Howes (2025) Multiobjective platform for autonomous property targeting and optimization of colloidal lead halide perovskite quantum dots. Chem. Mater. 37 (17), p. 6629–6641. External Links: ISBN 0897-4756, Document, Link Cited by: Introduction. [29] A. M. K. Nambiar, C. P. Breen, T. Hart, T. Kulesza, T. F. Jamison, and K. F. Jensen (2022) Bayesian optimization of computer-proposed multistep synthetic routes on an automated robotic flow platform. ACS Cent. Sci. 8 (6), p. 825–836. External Links: ISBN 2374-7943, Document, Link Cited by: Introduction. [30] H. Narayanan, J. A. Hinckley, R. Barry, B. Dang, L. A. Wolffe, A. Atari, Y. Tseng, and J. C. Love (2025) Accelerating cell culture media development using Bayesian optimization-based iterative experimental design. Nature Communications 16 (1), p. 6055. External Links: ISSN 2041-1723, Document, Link Cited by: Introduction. [31] N. Ndahiro, E. Ma, T. Bertalan, M. Donohue, Y. Kevrekidis, and M. Betenbaugh (2025) Integration of Bayesian optimization and solution thermodynamics to optimize media design for mammalian biomanufacturing. iScience 28 (8), p. 112944. External Links: ISSN 2589-0042, Document, Link Cited by: Introduction. [32] L. Papenmeier, N. Cheng, S. Becker, and L. Nardi (2025) Exploring Exploration in Bayesian Optimization. In The 41st Conference on Uncertainty in Artificial Intelligence, External Links: Link Cited by: Analysis of the optimization process, Analysis of the optimization process, Robustness analysis of hyperparameters and acquisition functions. [33] C. Qian, H. Xiong, and K. Xue (2021) Bayesian optimization using pseudo-points. In Proceedings of the Twenty-Ninth International Conference on International Joint Conferences on Artificial Intelligence, p. 3044–3050. Cited by: Introduction, Performance comparison with vanilla BO, Performance of PolyBO on synthetic benchmark functions. [34] J. P. Quek, A. A. Gaffoor, Y. X. Tan, T. R. M. Tan, Y. F. Chua, D. S. Z. Leong, A. S. Ali, and S. K. Ng (2024) Exploring cost reduction strategies for serum free media development. npj Science of Food 8 (1), p. 107. External Links: ISSN 2396-8370, Document, Link Cited by: Introduction, Introduction. [35] J. K. Salmon, M. A. Moraes, R. O. Dror, and D. E. Shaw (2011) Parallel random numbers: as easy as 1, 2, 3. In Proceedings of 2011 International Conference for High Performance Computing, Networking, Storage and Analysis, SC ’11, New York, NY, USA. External Links: ISBN 9781450307710, Link, Document Cited by: Numerical experiments using synthetic benchmark functions. [36] S. Skrivergaard, J. F. Young, N. Sahebekhtiari, C. Semper, M. Venkatesan, A. Savchenko, P. J. Stogios, M. Therkildsen, and M. K. Rasmussen (2023) A simple and robust serum-free media for the proliferation of muscle cells. Food Research International 172, p. 113194. External Links: ISSN 0963-9969, Document, Link Cited by: Introduction. [37] N. Srinivas, A. Krause, S. Kakade, and M. Seeger (2010) Gaussian process optimization in the bandit setting: no regret and experimental design. In Proceedings of the 27th International Conference on International Conference on Machine Learning, Vol. 1015–1022. Cited by: PolyBO algorithm. [38] Y. Taguchi, Y. Shibuya, Y. Hiki, T. Morikura, T. G. Yamada, and A. Funahashi (2026) Asynchronous batch Bayesian optimization with pipelining evaluations in experimental equipment-limited situations. SLAS Technology 37, p. 100396. External Links: ISSN 2472-6303, Document, Link Cited by: Performance comparison with vanilla BO. [39] G. Tom, S. P. Schmid, S. G. Baird, Y. Cao, K. Darvish, H. Hao, S. Lo, S. Pablo-García, E. M. Rajaonson, M. Skreta, N. Yoshikawa, S. Corapi, G. D. Akkoc, F. Strieth-Kalthoff, M. Seifrid, and A. Aspuru-Guzik (2024) Self-Driving Laboratories for Chemistry and Materials Science. Chemical Reviews 124 (16), p. 9633–9732. External Links: ISSN 0009-2665, Document, Link Cited by: Introduction. [40] A. Tripp, E. Daxberger, and J. M. Hernández-Lobato (2020) Sample-efficient optimization in the latent space of deep generative models via weighted retraining. In Advances in Neural Information Processing Systems, H. Larochelle, M. Ranzato, R. Hadsell, M.F. Balcan, and H. Lin (Eds.), Vol. 33, p. 11259–11272. External Links: Link Cited by: Introduction. [41] G. Walsh and E. Walsh (2022) Biopharmaceutical benchmarks 2022. Nature Biotechnology 40 (12), p. 1722–1760. External Links: ISSN 1546-1696, Document, Link Cited by: Introduction. [42] H. Wang, T. Fu, Y. Du, W. Gao, K. Huang, Z. Liu, P. Chandak, S. Liu, P. Van Katwyk, A. Deac, et al. (2023) Scientific discovery in the age of artificial intelligence. Nature 620 (7972), p. 47–60. Cited by: Introduction, Introduction. [43] J. Wang, H. Kwon, H. S. Kim, and B. Lee (2023) A neural network model for high entropy alloy design. npj Computational Materials 9 (1), p. 60. External Links: ISBN 2057-3960, Document, Link Cited by: Application to a real-world material composition optimization problem. [44] X. Wang, Y. Jin, S. Schmitt, and M. Olhofer (2023) Recent Advances in Bayesian Optimization. ACM Comput. Surv. 55 (13s). External Links: Document Cited by: Introduction. [45] Y. Xian, X. Ding, X. Jiang, Y. Zhou, J. Sun, D. Xue, and T. Lookman (2025) Unlocking the black box beyond Bayesian global optimization for materials design using reinforcement learning. npj Computational Materials 11 (1), p. 143. External Links: ISBN 2057-3960, Document, Link Cited by: Application to a real-world material composition optimization problem, Application to a real-world material composition optimization problem, Application to a real-world material composition optimization problem, Figure 7. [46] K. K. Yang, Z. Wu, and F. H. Arnold (2019) Machine-learning-guided directed evolution for protein engineering. Nature Methods 16 (8), p. 687–694. External Links: ISSN 1548-7105, Document, Link Cited by: Introduction. [47] Y. Yin, Y. Wang, and P. Li (2024) High-Dimensional Bayesian Optimization via Semi-Supervised Learning with Optimized Unlabeled Data Sampling. In International Conference on Machine Learning, p. 57085–57100. Cited by: Introduction, Performance comparison with vanilla BO, Performance of PolyBO on synthetic benchmark functions. [48] Z. Zanjani Foumani, M. Shishehbor, A. Yousefpour, and R. Bostanabad (2023) Multi-fidelity cost-aware Bayesian optimization. Computer Methods in Applied Mechanics and Engineering 407, p. 115937. External Links: ISSN 0045-7825, Document, Link Cited by: Discussion. [49] J. Zhao, R. Yang, S. Qiu, and Z. Wang (2024) Unleashing the Potential of Acquisition Functions in High-Dimensional Bayesian Optimization. Transactions on Machine Learning Research. Note: External Links: ISSN 2835-8856, Link Cited by: Robustness analysis of hyperparameters and acquisition functions. [50] H. Zheng, W. Xie, I. O. Ryzhov, and D. Xie (2023) Policy optimization in dynamic Bayesian network hybrid models of biomanufacturing processes. INFORMS Journal on Computing 35 (1), p. 66–82. Cited by: Introduction. Code availability The source code of this study is available from https://github.com/funalab/PolyBO. Acknowledgements The research was funded by JST CREST, Japan Grant Number JPMJCR21N1 to A.F. We are grateful to Asst. Prof. Yuki TSUKADA for helpful discussions. Author contributions H. S. : Data curation, Formal analysis, Investigation, Methodology, Software, Validation, Visualization, Writing–original draft, Writing–review & editing. Y. T. : Formal analysis, Methodology, Software, Visualization. K. M. : Formal analysis, Methodology, Software. Y. H. : Formal analysis, Writing–review & editing. T. M. : Conceptualization, Formal analysis, Methodology, Supervision, Visualization, Writing–review & editing. A. F. : Conceptualization, Formal analysis, Funding acquisition, Methodology, Project administration, Resources, Supervision, Visualization, Writing–review & editing. Competing Interests The authors declare no competing interests. Figure 1.: Overview of this study Schematic illustration of vanilla BO and PolyBO. (a) Vanilla BO trains a Gaussian process regression model using only the experimental data up to iteration t and calculates the acquisition function. Optimization is then performed by proposing one condition that maximizes the acquisition function, conducting an experiment under that condition, and obtaining the experimental data at iteration t+1t+1. (b) PolyBO combines vanilla BO with generation of regression-based pseudo-experimental data. The regression model is trained using the experimental data up to iteration t, and pseudo-experimental data are generated using the trained model. The pseudo-experimental data are pairs of explainable variables randomly sampled from the search space, and objective values are calculated by substituting the explainable variables into the regression model. The experimental data up to iteration t and the pseudo-experimental data are treated equally and combined into a dataset to train a Gaussian process regression model and calculate the acquisition function. Optimization is then performed as above. Figure 2.: Optimization processes using PolyBO and conventional methods The median (lines) and interquartile range (shading) of 30 runs of optimization from different initial conditions are shown for each method. D, dimensions; function IDs (F) are shown for 5 functions selected from the 24 functions. m′=10m =10, p=4p=4; TSBO was used only at D=20D=20. Table 1.: The number of iterations required for optimization by PolyBO and vanilla BO ID Benchmark function D=2D=2 D=5D=5 D=10D=10 D=20D=20 F01 Sphere (8) 17\,(8) – – (12) 49\,(12) F02 Separable ellipsoidal – – – (18) 39\,(18) F03 Rastrigin, original (58) 48\,(58) – (22) 27\,(22) (9) 23\,(9) F04 Büche–Rastrigin 65(–)65\,( ) 76(–)76\,( ) (38) 30\,(38) (19) 19\,(19) F05 Linear Slope (7) 8\,(7) (9) 16\,(9) (13) 44\,(13) 74(–)74\,( ) F06 Attractive sector 87(–)87\,( ) (29) 24\,(29) (33) 35\,(33) (16) 28\,(16) F07 Step ellipsoidal (52) 32\,(52) – 71(–)71\,( ) (14) 29\,(14) F08 Rosenbrock, original 59(18)59\,(18) – 73(–)73\,( ) (7) 26\,(7) F09 Rosenbrock, rotated 57(25)57\,(25) – (28) 27\,(28) (13) 22\,(13) F10 Ellipsoidal 91(–)91\,( ) – 68(–)68\,( ) (34) 39\,(34) F11 Discus 90(–)90\,( ) – (66) 41\,(66) (58) 45\,(58) F12 Bent cigar (34) 46\,(34) – – (14) 24\,(14) F13 Sharp ridge 67(45)67\,(45) – 65(–)65\,( ) 50(9)50\,(9) F14 Different Powers 62(37)62\,(37) – 72(–)72\,( ) (11) 30\,(11) F15 Rastrigin, rotated (36) 48\,(36) – (21) 28\,(21) (13) 13\,(13) F16 Weierstrass – 70(–)70\,( ) 64(–)64\,( ) 94(–)94\,( ) F17 Schaffer’s F7 59(48)59\,(48) – (20) 32\,(20) (21) 15\,(21) F18 Schaffer’s F7, moderately ill-conditioned – – (21) 26\,(21) (20) 16\,(20) F19 Composite Griewank–Rosenbrock function F8F2 (44) 31\,(44) – (37) 33\,(37) (9) 10\,(9) F20 Schwefel (46) 42\,(46) – (46) 31\,(46) (11) 30\,(11) F21 Gallagher’s Gaussian 101-me peaks 54(58)54\,(58) – (–) 36\,( ) 74(–)74\,( ) F22 Gallagher’s Gaussian 21-hi peaks 56(43)56\,(43) – 52(–)52\,( ) (53) 47\,(53) F23 Katsuura 58(75)58\,(75) – (36) 26\,(36) 63(–)63\,( ) F24 Lunacek bi-Rastrigin (62) 47\,(62) – (21) 26\,(21) (16) 12\,(16) Median 56.0 47.0 34.5 29.0 To quantify the reduction in optimization time, we calculated the number of iterations required for PolyBO to reach the simple regret reached by vanilla BO at 50 iterations. Each value is the median number of iterations at which PolyBO reached the reference simple regret, and the values in parentheses indicate the interquartile range. Values <50<50 are in bold, and values ≥50≥ 50 are in plain font. The maximum number of iterations was set to 100; – indicates that the value could not be calculated because the reference simple regret was not reached within 100 iterations. m′=10m =10, p=4p=4. Figure 3.: t-SNE visualization of optimization processes for F04 and F06 at D=20D=20 Circles indicate all sampled points used in the optimization, including the initial experimental points and the candidate points (t)x^(t) selected by each method at subsequent iterations. Color indicates the value of the objective function. A symlog scale was used for the color bar, and the threshold for the linear region was set to 10−310^-3 of the maximum absolute plotted value of the objective function. Red stars indicate the locations of the global optima. All results are shown for m′=10m =10 and p=4p=4. Figure 4.: Quantification of optimization behavior for F04 and F06 at D=20D=20 using OE and PE The median (lines) and interquartile range (shading) of 30 runs of optimization from different initial conditions are shown for each method. m′=10m =10, p=4p=4. Table 2.: Sensitivity analysis of the pseudo-experimental dataset size m′m in PolyBO ID m′=1m =1 m′=5m =5 m′=10m =10 m′=15m =15 m′=20m =20 m′=25m =25 m′=40m =40 m′=80m =80 m′=100m =100 m′=200m =200 F01 53(8)53\,(8) 52(14)52\,(14) (12) 49\,(12) 56(10)56\,(10) 50(10)50\,(10) 59(22)59\,(22) 57(25)57\,(25) 97(–)97\,( ) – – F02 (20) 43\,(20) (17) 46\,(17) (18) 39\,(18) (24) 41\,(24) (24) 42\,(24) (26) 42\,(26) (25) 42\,(25) 55(40)55\,(40) 58(46)58\,(46) 60(51)60\,(51) F03 (17) 31\,(17) (17) 29\,(17) (9) 23\,(9) (16) 25\,(16) (7) 23\,(7) (8) 25\,(8) (11) 29\,(11) (18) 28\,(18) (27) 34\,(27) (43) 42\,(43) F04 (16) 36\,(16) (20) 19\,(20) (19) 19\,(19) (18) 17\,(18) (18) 18\,(18) (17) 18\,(17) (10) 16\,(10) (14) 14\,(14) (15) 12\,(15) (14) 12\,(14) F05 65(40)65\,(40) 69(40)69\,(40) 74(–)74\,( ) – – – – – – – F06 (18) 39\,(18) (11) 28\,(11) (16) 28\,(16) (14) 31\,(14) (9) 27\,(9) (14) 27\,(14) (20) 31\,(20) (19) 35\,(19) (19) 39\,(19) (59) 47\,(59) F07 (27) 36\,(27) (12) 31\,(12) (14) 29\,(14) (13) 30\,(13) (13) 33\,(13) (13) 33\,(13) (17) 37\,(17) (26) 35\,(26) (19) 37\,(19) (17) 34\,(17) F08 (11) 39\,(11) (9) 34\,(9) (7) 26\,(7) (13) 25\,(13) (9) 25\,(9) (12) 24\,(12) (15) 28\,(15) (22) 27\,(22) (59) 45\,(59) 62(–)62\,( ) F09 (12) 37\,(12) (6) 25\,(6) (13) 22\,(13) (14) 19\,(14) (8) 15\,(8) (10) 14\,(10) (7) 10\,(7) (3) 8\,(3) (4) 8\,(4) (5) 5\,(5) F10 (23) 45\,(23) (42) 49\,(42) (34) 39\,(34) (23) 42\,(23) 59(42)59\,(42) (36) 46\,(36) 52(35)52\,(35) 73(–)73\,( ) 58(–)58\,( ) – F11 65(–)65\,( ) 56(78)56\,(78) (58) 45\,(58) 71(–)71\,( ) 68(–)68\,( ) 81(–)81\,( ) 76(–)76\,( ) 66(–)66\,( ) 94(–)94\,( ) 50(74)50\,(74) F12 (21) 37\,(21) (12) 26\,(12) (14) 24\,(14) (8) 23\,(8) (11) 18\,(11) (10) 19\,(10) (10) 15\,(10) (14) 19\,(14) (7) 19\,(7) (14) 18\,(14) F13 (7) 47\,(7) (13) 45\,(13) 50(9)50\,(9) 51(15)51\,(15) (15) 48\,(15) (10) 49\,(10) 54(22)54\,(22) 85(33)85\,(33) 91(–)91\,( ) – F14 (16) 39\,(16) (15) 29\,(15) (11) 30\,(11) (6) 29\,(6) (13) 32\,(13) (10) 33\,(10) (16) 31\,(16) (27) 35\,(27) (21) 42\,(21) 50(52)50\,(52) F15 (24) 24\,(24) (10) 12\,(10) (13) 13\,(13) (9) 13\,(9) (15) 13\,(15) (11) 13\,(11) (10) 12\,(10) (15) 17\,(15) (14) 15\,(14) (10) 14\,(10) F16 50(62)50\,(62) 67(–)67\,( ) 94(–)94\,( ) 55(–)55\,( ) – – – – – – F17 (22) 26\,(22) (20) 18\,(20) (21) 15\,(21) (23) 15\,(23) (23) 13\,(23) (11) 14\,(11) (14) 12\,(14) (12) 10\,(12) (9) 14\,(9) (8) 11\,(8) F18 (29) 26\,(29) (23) 21\,(23) (20) 16\,(20) (28) 15\,(28) (13) 18\,(13) (17) 13\,(17) (13) 13\,(13) (14) 15\,(14) (12) 14\,(12) (16) 12\,(16) F19 (14) 15\,(14) (10) 10\,(10) (9) 10\,(9) (7) 9\,(7) (12) 10\,(12) (6) 8\,(6) (8) 9\,(8) (9) 9\,(9) (8) 8\,(8) (10) 6\,(10) F20 (11) 43\,(11) (10) 31\,(10) (11) 30\,(11) (9) 28\,(9) (11) 29\,(11) (9) 29\,(9) (14) 35\,(14) (35) 45\,(35) 63(–)63\,( ) – F21 (–) 48\,( ) (–) 38\,( ) 74(–)74\,( ) 82(–)82\,( ) 94(–)94\,( ) 56(–)56\,( ) 85(–)85\,( ) – 74(–)74\,( ) – F22 69(–)69\,( ) 63(–)63\,( ) (53) 47\,(53) 52(72)52\,(72) 56(54)56\,(54) 54(46)54\,(46) (46) 42\,(46) 65(48)65\,(48) 58(32)58\,(32) 69(24)69\,(24) F23 (42) 38\,(42) 70(–)70\,( ) 63(–)63\,( ) 61(–)61\,( ) 60(–)60\,( ) (61) 36\,(61) (–) 44\,( ) (51) 35\,(51) 51(–)51\,( ) (–) 41\,( ) F24 (36) 31\,(36) (25) 16\,(25) (16) 12\,(16) (12) 14\,(12) (7) 10\,(7) (15) 11\,(15) (15) 12\,(15) (12) 10\,(12) (13) 9\,(13) (8) 5\,(8) Median 38.5 30.8 29.0 29.0 28.0 27.8 30.8 34.5 38.5 34.0 We calculated the number of iterations required for PolyBO to reach the simple regret reached by vanilla BO at 50 iterations. Each value is the median number of iterations at which PolyBO reached the reference simple regret, and the values in parentheses indicate the interquartile range. Values <50<50 are in bold, and values ≥50≥ 50 are in plain font. The maximum number of iterations was set to 100; – indicates that the value could not be calculated because the reference simple regret was not reached within 100 iterations. D=20D=20, p=4p=4. Table 3.: Sensitivity analysis of the maximum polynomial degree p in PolyBO ID p=1p=1 p=2p=2 p=3p=3 p=4p=4 p=5p=5 p=6p=6 p=7p=7 p=8p=8 F01 72(9)72\,(9) (6) 48\,(6) – (12) 49\,(12) (6) 48\,(6) 63(36)63\,(36) – – F02 55(22)55\,(22) (19) 43\,(19) 67(39)67\,(39) (18) 39\,(18) (19) 43\,(19) (19) 39\,(19) (33) 48\,(33) 50(28)50\,(28) F03 (8) 38\,(8) (10) 30\,(10) (16) 30\,(16) (9) 23\,(9) (10) 30\,(10) (18) 21\,(18) (14) 22\,(14) (15) 27\,(15) F04 (46) 48\,(46) (19) 26\,(19) (36) 36\,(36) (19) 19\,(19) (19) 26\,(19) (22) 19\,(22) (23) 24\,(23) (18) 19\,(18) F05 55(7)55\,(7) 71(25)71\,(25) 73(–)73\,( ) 74(–)74\,( ) 71(25)71\,(25) 78(–)78\,( ) 75(–)75\,( ) 79(–)79\,( ) F06 (5) 43\,(5) (15) 38\,(15) (14) 32\,(14) (16) 28\,(16) (15) 38\,(15) (16) 30\,(16) (21) 30\,(21) (16) 31\,(16) F07 (11) 45\,(11) (18) 37\,(18) (25) 32\,(25) (14) 29\,(14) (18) 37\,(18) (18) 32\,(18) (13) 32\,(13) (16) 32\,(16) F08 52(17)52\,(17) (9) 32\,(9) (19) 45\,(19) (7) 26\,(7) (9) 32\,(9) (11) 28\,(11) (36) 41\,(36) (20) 29\,(20) F09 – (7) 26\,(7) (18) 38\,(18) (13) 22\,(13) (7) 26\,(7) (11) 22\,(11) (16) 26\,(16) (13) 17\,(13) F10 83(–)83\,( ) 53(39)53\,(39) 100(–)100\,( ) (34) 39\,(34) 53(39)53\,(39) (23) 38\,(23) 57(–)57\,( ) (46) 44\,(46) F11 84(–)84\,( ) 73(–)73\,( ) 74(–)74\,( ) (58) 45\,(58) 73(–)73\,( ) 52(–)52\,( ) 84(–)84\,( ) – F12 (12) 38\,(12) (12) 28\,(12) (13) 24\,(13) (14) 24\,(14) (12) 28\,(12) (12) 24\,(12) (13) 23\,(13) (10) 24\,(10) F13 70(8)70\,(8) (8) 47\,(8) 72(26)72\,(26) 50(9)50\,(9) (8) 47\,(8) 52(14)52\,(14) – 82(–)82\,( ) F14 53(21)53\,(21) (15) 34\,(15) (26) 39\,(26) (11) 30\,(11) (15) 34\,(15) (13) 31\,(13) (21) 39\,(21) (25) 36\,(25) F15 (13) 32\,(13) (15) 13\,(15) (12) 17\,(12) (13) 13\,(13) (15) 13\,(15) (12) 10\,(12) (11) 11\,(11) (5) 9\,(5) F16 83(–)83\,( ) – – 94(–)94\,( ) 91(–)91\,( ) (59) 47\,(59) 93(–)93\,( ) 63(61)63\,(61) F17 (19) 23\,(19) (17) 21\,(17) (17) 14\,(17) (21) 15\,(21) (14) 17\,(14) (16) 19\,(16) (16) 10\,(16) (15) 14\,(15) F18 (16) 23\,(16) (22) 21\,(22) (18) 10\,(18) (20) 16\,(20) (14) 13\,(14) (17) 15\,(17) (20) 11\,(20) (31) 20\,(31) F19 (9) 42\,(9) (22) 28\,(22) (13) 16\,(13) (9) 10\,(9) (8) 10\,(8) (11) 10\,(11) (7) 9\,(7) (7) 8\,(7) F20 61(17)61\,(17) (14) 35\,(14) (16) 41\,(16) (11) 30\,(11) (15) 39\,(15) (12) 28\,(12) (31) 47\,(31) (18) 37\,(18) F21 (–) 41\,( ) (18) 32\,(18) 53(–)53\,( ) 74(–)74\,( ) (–) 40\,( ) (–) 44\,( ) (–) 41\,( ) (69) 44\,(69) F22 52(58)52\,(58) 51(36)51\,(36) 77(–)77\,( ) (53) 47\,(53) (48) 48\,(48) 51(75)51\,(75) (42) 31\,(42) (41) 38\,(41) F23 70(–)70\,( ) (84) 35\,(84) 52(–)52\,( ) 63(–)63\,( ) (86) 38\,(86) 58(–)58\,( ) (76) 46\,(76) (54) 29\,(54) F24 (30) 37\,(30) (26) 24\,(26) (27) 40\,(27) (16) 12\,(16) (11) 13\,(11) (16) 12\,(16) (14) 14\,(14) (10) 10\,(10) Median 51.5 33.5 39.3 29.0 37.3 30.3 31.3 30.0 We calculated the number of iterations required for PolyBO to reach the simple regret reached by vanilla BO at 50 iterations. Each value is the median number of iterations at which PolyBO reached the reference simple regret, and the values in parentheses indicate the interquartile range. Values <50<50 are in bold and values ≥50≥ 50 are in plain font. The maximum number of iterations was set to 100; – indicates that the value could not be calculated because the reference simple regret was not reached within 100 iterations. D=20D=20, m′=10m =10. Figure 5.: Sensitivity analysis of the acquisition function in PolyBO EI or GP-UCB was used as the acquisition function. The median (lines) and interquartile range (shading) of 30 runs of optimization from different initial conditions are shown for each method. D=20D=20, m′=10m =10, and p=4p=4. Figure 6.: Ablation study of the update mechanism The update mechanism was ablated by either not discarding previously generated pseudo-experimental data (w/o reset) or increasing the pseudo-experimental dataset size (m′m -scaling). The median (lines) and interquartile range (shading) of 30 runs of optimization from different initial conditions are shown for each PolyBO variant. D=20D=20, m′=10m =10, and p=4p=4. Figure 7.: Real-world material composition optimization with PolyBO and two conventional methods The median (lines) and interquartile range (shading) are shown. The results of 30 (PolyBO and BO-EI) or 24 (RL-DQN) independent runs with different initial conditions are shown. The initial conditions in PolyBO were the same 20 conditions as used in BO-EI by Xian et al. [45]. The results of BO-EI and RL-DQN were obtained by reproducing the published code [45]. D=10D=10, m′=10m =10, and p=4p=4. Supplementary Figure 1.: Optimization processes for 24 benchmark functions at D=2D=2 The median (lines) and interquartile range (shading) of 30 runs of optimization from different initial conditions are shown for each method. Problem setting: m′=10m =10, p=4p=4. Supplementary Figure 2.: Optimization processes for 24 benchmark functions at D=5D=5 The median (lines) and interquartile range (shading) of 30 runs of optimization from different initial conditions are shown for each method. Problem setting: m′=10m =10, p=4p=4. Supplementary Figure 3.: Optimization processes for 24 benchmark functions at D=10D=10 The median (lines) and interquartile range (shading) of 30 runs of optimization from different initial conditions are shown for each method. Problem setting: m′=10m =10, p=4p=4. Supplementary Figure 4.: Optimization processes for 24 benchmark functions at D=20D=20 The median (lines) and interquartile range (shading) of 30 runs of optimization from different initial conditions are shown for each method. Problem setting: m′=10m =10, p=4p=4. Supplementary Figure 5.: UMAP visualization of optimization processes for F04 and F06 at D=20D=20 Circles indicate all sampled points used in the optimization, including the initial experimental points and the candidate points (t)x^(t) selected by each method at subsequent iterations. Color indicates the value of the objective function. A symlog scale was used for the color bar, and the threshold for the linear region was set to 10−310^-3 of the maximum value of the absolute plotted objective function. Red stars indicate the locations of the global optima. m′=10m =10, p=4p=4. Supplementary Figure 6.: t-SNE visualization of optimization processes for F04 and F06 at D=2D=2 Circles indicate all sampled points used in the optimization, including the initial experimental points and the candidate points (t)x^(t) selected by each method at subsequent iterations. Color indicates the value of the objective function. A symlog scale was used for the color bar, and the threshold for the linear region was set to 10−310^-3 of the maximum absolute plotted objective function value. Red stars indicate the locations of the global optima. m′=10m =10, p=4p=4. Supplementary Figure 7.: t-SNE visualization of optimization processes for F04 and F06 at D=5D=5 Circles indicate all sampled points used in the optimization, including the initial experimental points and the candidate points (t)x^(t) selected by each method at subsequent iterations. Color indicates the value of the objective function. A symlog scale was used for the color bar, and the threshold for the linear region was set to 10−310^-3 of the maximum absolute plotted objective function value. Red stars indicate the locations of the global optima. m′=10m =10, p=4p=4. Supplementary Figure 8.: Entire search space after dimensionality reduction using t-SNE A symlog scale was used for the color bar, and the threshold for the linear region was set to 10−310^-3 of the value of the maximum absolute plotted objective function. Color indicates the value of the objective function. Red stars indicate the locations of the global optima. m′=10m =10, p=4p=4. Supplementary Figure 9.: Quantification of optimization behavior for F04 and F06 at D=2D=2 using OE and PE The median (lines) and interquartile range (shading) of 30 runs of optimization from different initial conditions are shown for each method. m′=10m =10, p=4p=4. Supplementary Figure 10.: Quantification of optimization behavior for F04 and F06 at D=5D=5 using OE and PE The median (lines) and interquartile range (shading) of 30 runs of optimization from different initial conditions are shown for each method. m′=10m =10, p=4p=4. Supplementary Algorithm 1 PolyBO algorithm for a discrete constrained search space 0: Discrete grid interval δ, bounds [li,ui](i=1,…,d)[l_i,u_i]\ (i=1,…,d) for each variable, equality constraint ∑i=1dxi=1 _i=1^dx_i=1 1: Uniformly sample k explainable variables from the discrete grid set δG_δ satisfying ∑i=1dxi=1 _i=1^dx_i=1 and li≤xi≤uil_i≤ x_i≤ u_i, and construct the initial experimental dataset 0:0=(i(0),yi(0))i=1kD_0:0=\(x^(0)_i,y^(0)_i)\_i=1^k. 2: for t=1,2,…t=1,2,… do 3: Train pthp^th degree polynomial regression model fp()f_p(x) using dataset 0:t−1D_0:t-1. 4: Uniformly sample m′m explainable variables ′j(t)j=1m′\x ^(t)_j\_j=1^m from the discrete grid set δG_δ. 5: Compute yj′(t)=fp(′j(t))y (t)_j=f_p(x ^(t)_j) for j=1,…,m′j=1,…,m , and construct pseudo-experimental data t′=(′j(t),yj′(t))j=1m′D_t =\(x ^(t)_j,y (t)_j)\_j=1^m . 6: Update the Gaussian process regression model with 0:t−1∪t′D_0:t-1 _t . 7: Compute the acquisition function α()α(x) using the updated Gaussian process regression model. 8: Optimize α()α(x) in a continuous space under the constraint ∑i=1dxi=1 _i=1^dx_i=1 using the SLSQP method from multiple initial conditions. 9: For each continuous optimum, enumerate adjacent lower and upper discrete grid points in each dimension, and extract only combinations satisfying ∑i=1dxi=1 _i=1^dx_i=1 as discrete candidates. 10: Re-evaluate α()α(x) for all discrete candidates and select the unobserved candidate with the maximum value as (t)x^(t). 11: Conduct an experiment at (t)x^(t) to obtain y(t)y^(t). 12: Update the dataset 0:t=0:t−1∪((t),y(t))D_0:t=D_0:t-1∪\(x^(t),y^(t))\. 13: end for Supplementary Figure 11.: Peak memory usage and number of parameters at different maximum polynomial degrees Peak memory usage (a) was measured to evaluate memory complexity.