Data-Driven Prior Construction in Hilbert Spaces for Bayesian Optimization

1. Introduction

Optimization of expensive-to-evaluate black-box functions is a problem found in a wide range of technological and engineering applications [1,2,3]. Two broad application fields include the tuning of Machine Learning (ML) and general Artificial Intelligence (AI) models[4,5,6], as well as multiple forms of engineering design optimization, such as shape, topology, and material optimization, spanning from structural to electronic engineering [7,8,9]. In both fields, complex computational models are required to compute the objective or performance function, resulting in elevated computation time for each evaluation of the objective. Additionally, these models act as black boxes, meaning that general properties of the underlying model (such as its derivative or gradient) are not available.

Bayesian Optimization (BO) is an optimization framework particularly well-suited to address the aforementioned black-box, expensive-to-evaluate optimization problem [10,11,12,13,14]. It harnesses the efficacy of a surrogate model to approximate the objective function. It utilizes a Bayesian predictor (an Acquisition function) of promising enrichment points that are likely to improve the prediction of the surrogate in regions of interest, that is, regions where the optimum is likely to be found based on the limited information available.

Ideally, an instance of BO balances exploration and exploitation [15] by combining the prediction and uncertainty of the surrogate, usually a Gaussian Process ($\mathcal{GP}$), with the acquisition function [16,17].

Since BO is Bayesian in nature, its performance depends on the quality of the prior structure used in the model. No single prior can be ideal to capture every type of objective function (NFL principle), and a misguided prior may lead to unfeasible requirements in terms of functional evaluations, or to failure to converge towards a feasible optimum. In this sense, several versions of BO have explored ways to incorporate adapted priors into the basic BO algorithm [18,19], to accelerate convergence. In particular, Bayesian optimization with a prior for the optimum (BOPrO) [20] and related methods introduce a belief based on expert knowledge, which allows certain regions of the input space to be favoured as more promising. Strictly speaking, these approaches codify a form of prior that stems from the user’s previous experience and guide the acquisition function to explore such regions.

This article presents an algorithm that combines the advantages of classical Bayesian optimisation with prior distributions constructed from an initial experimental design (DoE, Design of Experiment), based on Uncertainty Quantification (UQ) techniques. While both the metamodel and prior distributions utilize the same DoE, they encode complementary information about the problem. While the metamodel seeks to approximate the objective function across the entire domain, taking advantage of spatial correlations between samples, the prior distributions quantify the search space in a Bayesian sense by identifying regions where the objective is likely to attain good or bad performance contingent on the observations.

In engineering contexts, expert priors are often unavailable or difficult to justify, particularly in early design phases. In such cases, investing in a small exploratory DoE enables data-driven prior construction that can guide the optimization process efficiently from the outset, as demonstrated in this work.

The organization of this article is the following: Section 2 reviews the fundamentals of Bayesian optimization and the original BOPrO scheme. Section 3 details the construction of the UQ-based performance prior and its integration into the complete algorithm. Section 4 presents numerical results from six academic test functions and a case study from an engineering application, highlighting the applicability of the proposed method. Finally, Section 5 summarizes the conclusions and outlines future research directions.

2. Background

We are interested in solving the following global optimization problem:

$${\mathbf{x}}^{*}=argmin\left\{f\left(\mathbf{x}\right):\mathbf{x}\in \mathcal{X}\subset {\mathbb{R}}^d\right\},$$

(1)

where $f:{\mathbb{R}}^d\to \mathbb{R}$ is the function to be minimized.

Sequential Bayesian estimation gives rise to a global optimization framework, commonly referred to as BO, which is particularly suited for problems such as the one described above. BO relies on the fundamental principle of Bayesian updating. A typical iteration of a Bayesian Optimization Algorithm (BOA) proceeds as follows [21]:

• Generate a surrogate model of the objective function f, usually under the form of a Stochastic Process,
• Select a new point to evaluate the true objective function,
• Update the surrogate model using the Bayesian Update, the pair of the new point and the value of the objective function in it.

2.1. SBO

The most commonly used surrogate model in BOA is the $\mathcal{GP}$[14,15,16,17,22,23]. A $\mathcal{GP}$ generalizes the multivariate Gaussian distribution to an infinite-dimensional stochastic process, where any finite collection of function evaluations follows a joint multivariate normal distribution. Just as a standard Gaussian distribution describes a random variable by its mean and covariance, a $\mathcal{GP}$ defines a distribution over functions, fully specified by a mean function $\mathbf{m}\left(\mathbf{x}\right)$ and a covariance (or kernel) function $\mathbf{K}(\mathbf{x},{\mathbf{x}}^{\prime })$ [15]:

$$\widehat{f}\sim \mathcal{G}\mathcal{P}\left(\mathbf{m},\mathbf{K}\right),$$

(2)

where the mean function is often set to zero, $\mathbf{m}\left(\mathbf{x}\right)=0$, and the covariance matrix $\mathbf{K}$ has elements $K_{i,j}:=k({\mathbf{x}}_i,{\mathbf{x}}_j)$, with k denoting a positive semi-definite kernel function. Equation (2) represents the prior distribution $p\left(\widehat{f}\right)$ induced by the $\mathcal{GP}$, encoding assumptions about the smoothness and structure of the underlying objective function.

In terms of Bayesian update, the choice of a $\mathcal{GP}$ greatly simplifies the procedure, as the $\mathcal{GP}$ is conjugate to itself. This property ensures that the posterior distribution remains a $\mathcal{GP}$, and that the update merely involves a straightforward modification of the mean and covariance functions [21]. More precisely, we assume that a set of observations, ${\mathcal{D}}_t={\left\{\left({\mathbf{x}}_i,y_i\right)\right\}}_{i=1}^t$, has already been collected from previous iterations, where the outputs are noise-free evaluations of a deterministic function, i.e., $y_i=f\left({\mathbf{x}}_i\right)$. Let, $\mathbf{y}={\left[y_1,\dots ,y_t\right]}^T$, denote the vector of corresponding observations. Then, Bayes’ theorem yields the posterior distribution [15]:

$$p\left(\widehat{f}\mid {\mathcal{D}}_t\right)\propto p\left(\mathbf{y}\mid \widehat{f}\right)·p\left(\widehat{f}\right),$$

(3)

the posterior distribution $\widehat{f}\mid {\mathcal{D}}_t\sim \mathcal{GP}\left({\mu }_t,{\sigma }_t^2\right)$, which remains a $\mathcal{GP}$ due to the aforementioned conjugacy. This updated distribution leads to the following expression for the predictive distribution at any next evaluation point $\mathbf{x}$, which follows a univariate Gaussian distribution [15]:

$$\widehat{f}\sim \mathcal{N}\left(\mu ,{\sigma }^2\right),$$

(4)

with the predictive mean and variance given by:

$$\begin{array}{cc}\hfill \mu \left(\mathbf{x}\right)& ={\mathbf{k}}^T{\mathbf{K}}^{-1}\mathbf{y},\hfill \\ \hfill {\sigma }^2\left(\mathbf{x}\right)& =k(\mathbf{x},\mathbf{x})-{\mathbf{k}}^T{\mathbf{K}}^{-1}\mathbf{k},\hfill \end{array}$$

(5)

where, $\mathbf{k}=\left[k(\mathbf{x},{\mathbf{x}}_1),\dots ,k(\mathbf{x},{\mathbf{x}}_t)\right]$. To select the next point $\mathbf{x}$, in BOA one must solve an auxiliary optimization problem involving the maximization of the acquisition function $\alpha$. Among the popular acquisition functions, we can cite the Expected Improvement (EI). Let $f_{★}$ denote the best (i.e., lowest) function value observed up to iteration t, and assume the predictive distribution at $\mathbf{x}$ follows Equations (4) and (5). Then, the EI at a candidate point $\mathbf{x}$ is defined as follows [21]:

$${\alpha }_{EI}\left(\mathbf{x}\right)=\mathbb{E}\left[u_{EI}\left(\mathbf{x}\right)\mid {\mathcal{D}}_t\right],\mathrm{with}{u}_{EI}\left(\mathbf{x}\right)=max\left\{0,f_{★}-\widehat{f}\right\}.$$

(6)

A second popular acquisition function is the Probability of Improvement (PI):

$${\alpha }_{PI}\left(\mathbf{x}\right)=\mathbb{P}\left(\widehat{f}<f_{★}\mid {\mathcal{D}}_t\right),$$

(7)

Notice that:

$${\alpha }_{PI}\left(\mathbf{x}\right)=\mathbb{E}\left[u_{PI}\left(\mathbf{x}\right)\mid {\mathcal{D}}_t\right],\mathrm{with}{u}_{PI}\left(\mathbf{x}\right)=\left\{\begin{array}{cc}1\hfill & \mathrm{if}\widehat{f}<f_{★}\hfill \\ 0\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(8)

Another advantage of $\mathcal{GP}$ is that they allow the obtaining of closed analytical expressions for the EI and PI. Let $\mathsf{\Phi }$ and $\phi$, respectively, be the cumulative distribution function and the probability density function of the standard normal distribution $\mathcal{N}(0,1)$. Since the predictive distribution at a candidate point $\mathbf{x}$ is Gaussian with mean $\mu$ and standard deviation $\sigma$, the EI can be calculated as follows [21]:

$${\alpha }_{EI}\left(\mathbf{x}\right)=\sigma \left(\mathbf{x}\right)\left(z\left(\mathbf{x}\right)\mathsf{\Phi }\left(\right(\mathbf{x}\left)\right)+\phi \left(\right(\mathbf{x})\right),$$

(9)

with,

$$z\left(\mathbf{x}\right)={\displaystyle \frac{f_{★}-\mu \left(\mathbf{x}\right)}{\sigma \left(\mathbf{x}\right)}}.$$

(10)

Similarly, the PI can also be expressed as:

$${\alpha }_{PI}\left(\mathbf{x}\right)=\mathsf{\Phi }\left(z\left(\mathbf{x}\right)\right).$$

(11)

Throughout this article, we will refer to Standard Bayesian Optimization (SBO) as the Bayesian optimization framework where the surrogate model is a $\mathcal{GP}$ and the acquisition function is the EI.

To illustrate how the SBO works in practice, we will consider a one-dimensional synthetic test function proposed by Gramacy and Lee [24]. The objective is to find the global minimum of this function:

Figure 1. Illustration of the SBO optimization process. Figure (a) shows the shape of the function we want to optimize, and Figure (b) shows the mean (prediction) of the surrogate model with the associated uncertainty, as well as the structure of the EI and its maximum.

Figure 1. Illustration of the SBO optimization process. Figure (a) shows the shape of the function we want to optimize, and Figure (b) shows the mean (prediction) of the surrogate model with the associated uncertainty, as well as the structure of the EI and its maximum.

2.2. BOPrO

Bayesian optimization with a prior for the optimum (BOPrO), proposed by Souza et al. [20], is a recently developed version of BO based on the tree-structured Parzen estimator (TPE) [25], which allows users to introduce intuitive domain knowledge into the optimization process. These priors do not modify the surrogate model (e.g., a $\mathcal{GP}$). However, they are combined with the probability of improvement (PI) to form what is called a pseudo-posterior, more precisely, a posterior pseudo-density because it is not necessarily normalized (does not necessarily integrate to 1 in the prescribed domain), and acts as a weighting function inspired by Bayesian updating. This function is used in the acquisition function, in particular in the calculation of the modified expected improvement, which implies that the prior directly rebalances the weights assigned to the different regions of the search space in the EI evaluation. Thus, regions favoured by the prior are more likely to be explored at the beginning of the iteration process. As the iterations progress, the weight of the prior decreases, and the influence of the predictive model dominates the search.

The mathematical expression corresponding to this modified version of the expected improvement criterion, adopted by BOPrO to effectively integrate the pseudo-posterior into the process of selecting new points, is formally derived from the TPE framework. It reads:

$$E{I}_{f_{\gamma }}\left(\mathbf{x}\right)\propto {\left(\gamma +{\displaystyle \frac{b\left(\mathbf{x}\right)}{g\left(\mathbf{x}\right)}}(1-\gamma )\right)}^{-1},$$

(12)

where g and b are non-normalized density functions describing, respectively, the promising and less promising areas of the search space. The $\gamma$ value acts as a quantile threshold that determines at what point a function value can be considered “good”.

The main contribution of BOPrO is the practical implementation of this performance pseudo-posterior. In the case of $g\left(\mathbf{x}\right)$, a prior distribution $p_g\left(\mathbf{x}\right)$ is integrated, reflecting the user’s belief about promising regions, together with a probability derived from the model ${\mathcal{M}}_g\left(\mathbf{x}\right)$, which represents the possibility, under the surrogate $\mathcal{GP}$-that the value of the function $\mathbf{x}$ falls below a threshold $f_{\gamma }$. Which we can express as follows:

$$g\left(\mathbf{x}\right)\propto p_g\left(\mathbf{x}\right)·{\mathcal{M}}_g{\left(\mathbf{x}\right)}^{{\displaystyle \frac{t}{\beta }}},$$

(13)

The density for the unpromising regions is given by:

$$b\left(\mathbf{x}\right)\propto p_b\left(\mathbf{x}\right)·{\mathcal{M}}_b{\left(\mathbf{x}\right)}^{{\displaystyle \frac{t}{\beta }}},$$

(14)

where $p_b\left(\mathbf{x}\right)=1-p_g\left(\mathbf{x}\right)$, represents the prior for bad regions, t is the current iteration and $\beta$ hyperparameter defines the balance between prior and model, with higher $\beta$ values giving more importance to the prior and requiring more data to overrule it [20].

The probability ${\mathcal{M}}_g$ is defined as follows:

$${\mathcal{M}}_g\left(\mathbf{x}\right)=\mathsf{\Phi }\left({\displaystyle \frac{f_{\gamma }-\mu \left(x\right)}{\sigma \left(x\right)}}\right),$$

(15)

and the probability for unpromising regions is given by:

$${\mathcal{M}}_b\left(\mathbf{x}\right)=1-{\mathcal{M}}_g\left(\mathbf{x}\right).$$

(16)

Algorithm 1 describes in detail the optimization process carried out by BOPrO.

Algorithm 1: BOPrO Algorithm. ${\mathcal{D}}_t$ keeps track of all function evaluations so far: ${(x_i,y_i)}_{i=1}^t$ [20].

3. Data-Driven Prior Construction in Hilbert Spaces for Bayesian Optimization

3.1. Preliminaries

The BOPrO framework is a technical development adapted for machine learning applications, with the properties of the algorithm crafted explicitly for this type of problem. The performance priors are established in a way that facilitates use by an ML expert with domain knowledge in hyperparameter tuning. The goal is to provide a practical tool to codify this expert knowledge into a prior that: 1) describes the belief that certain regions of the search space are more (respectively, less) promising; 2) can be combined with the information available, both from model evaluations and from surrogate model predictions. This formulation poses significant limitations when the domain of application is different, such as structural design or general engineering design optimization. Notably, in such applications, expert knowledge may not be available, or sufficiently reliable, or the expert knowledge may be applied at a different stage of the process.

Our proposed method takes some of the central ideas in BOPrO (and the TPE algorithm). In particular, we retain the adapted version of the EI acquisition function and the structure of the priors over the search space, informing the early stages of the optimization iteration process. However, unlike the previous work, we replace the expert knowledge codified in the performance priors by systematic and reproducible distributions, learned directly from initial observations. To achieve this, we adopt the approach proposed by Souza and Fabro [26], which provides distributions that can be considered priors in a Bayesian procedure, particularly effective when data is limited.

3.2. Construction of Priors

The Hilbert Approach (HA), which can be viewed as an extension of Polynomial Chaos Approximations (PCA), can be applied under two basic conditions: is assumed to satisfy some basic conditions:

• The uncertainties can be modeled by a random vector U and the variable of interest is $X=X\left(U\right)$ ;
• the random variable U takes its values on an interval $\mathsf{\Omega }$ of real numbers (or, more generally, on a bounded subset $\mathsf{\Omega }$ of ${\mathbb{R}}^n$);
• X belongs to the space of square summable random variables: $E\left(X^2\right)<\infty$, id est, X has finite moment of order 2;
• some statistical information about the couple $(U,X)$ is available.

These assumptions ensure that X can be represented by an infinite expansion involving a convenient orthonormal Hilbert basis ${\left\{{\phi }_i\right\}}_{i\in \mathbb{N}}\subset L^2\left(\mathsf{\Omega }\right)$, allowing one to represent X as [26,27]:

$$X=\sum _{i\in \mathbb{N}}{x}_i{\phi }_i\left(U\right),$$

(17)

As previously observed, U is a suitable random variable used to construct a stochastic representation of X. Nevertheless, in practical situations, U can be an artificial variable, conveniently chosen [21,28]. Artificial variables are used, namely, when not any information is available about the source of randomness: then, U is chosen by the user and it is mandatory to create an artificial correlation between samples of the artificial variable and samples of X. For instance,we can range both in an increasing order: $X_1\le X_2\le \dots \le X_n$, and $U_1\le$$U_2\le \dots \le U_n$ - this generates a non-negative covariance between X and the artificial U.

In practice, the expansion is truncated at order k, yielding an approximation :

$$X=\sum _{i\in \mathbb{N}}{x}_i{\phi }_i\left(U\right)\approx P_{k}X=\sum _{i=0}^k{x}_i{\phi }_i\left(U\right),$$

(18)

The finite representation $P_{k}X$ is determined by finding the coefficients $\mathbf{x}={\left(x_0,x_1,\dots ,x_k\right)}^t$. A popular method for the determination of $\mathbf{x}$ is collocation, which considers a sample $\mathcal{U}=\left(U_1,\dots ,U_n\right)$ and solves the linear system [28]:

$$\mathit{M}\mathit{x}=\mathit{N},{\mathit{M}}_{ij}={\phi }_j\left(U_i\right),N_i=X_i,1\le i\le n,0\le j\le k$$

(19)

The reader can find in the literature works using such an approach. In general, $k+1<n$ (the system is overdetermined) and Equation (19) must be solved by a least squares approach. Once $P_{k}X$ is determined, the representation can be used to generate the distribution of X. For instance, we can generate a large sample of variates from U (eventually, a deterministic one),

$${\mathcal{U}}_g=\left(U_1,\dots ,U_{ng}\right)$$

(20)

and use it to generate a large sample ${\mathcal{X}}_g$ from X as $X_i=P_{k}X\left(U_i\right)$, $i=1,\dots ,ng$. Then, ${\chi }_g$ is used to estimate the cumulative distribution function (CDF) and probability density function (PDF) of X, by using the empirical CDF of the large sample.

From the Bayesian standpoint, estimates can be obtained by minimizing a loss function defined over the estimated PDF. This approach can be used to define a univariate parametric prior distribution $p_{\widehat{\theta }}\left(x\right)$, where the parametric family (e.g., Gaussian, log-normal) is selected by the user.

In the case of a multidimensional input, we assume independence between dimensions, so that the a priori multivariate distribution is constructed as the product of the univariate margins:

$$p\left(x\right)=\prod _{j=1}^d{p}_{{\widehat{\theta }}_j}\left(x_j\right).$$

(21)

This independence hypothesis, although simplistic and not necessarily valid in all cases, is adopted here as a reasonable working assumption to facilitate model construction and advance the validation of the proposed ideas.

The full optimization procedure is detailed in Algorithm (2).

Let’s return to the practical case defined in Section 2, this time applying our proposed method. The objective is to visualise the interaction between the different elements of the proposed method: the surrogate model, the acquisition function, the prior distributions, and the update points. Using the same DoE, we analyse how the behaviour of the acquisition function differs between our approach and SBO. The prior distributions we propose allow us to identify regions of the search space that, although considered promising by the metamodel, receive a low probability (see Figure 2(b)). As a result, the algorithm avoids exploring these regions and selects a new point different from the one that the traditional approach would have chosen (see Figure 2(f)).

Algorithm 2: Data-Driven Prior Construction in Hilbert Spaces for Bayesian Optimization

4. Results

In this section, we test our proposed method on six test functions and one engineering application. Our results are contrasted with SBO as a baseline. The test functions have been chosen to explore the general behavior of the method, even though they do not constitute expensive black-box functions. The application is a good representative of the kind of problem where BO is expected to excel.

Our focus here is exploratory, as we seek to identify interesting behavior of the proposed method against SBO, and with varied objective functions. We have thus selected performance metrics that reflect a mix of robustness(success rate), reliability (proximity to the optimum), and convergence speed (simple regret plots); moreover, particular attention is given to the variability of the method given different instances of the DoE, leading us to report extremes, averages and variances of each instance of our method for several runs.

4.1. Test Design and Technical Aspects

Our approach is evaluated by identifying the global minima of a set of test functions from the Virtual Library of Simulation Experiments [24].

Table A1 in Appendix A shows the mathematical expression, the problem size, the search space, the known global minimum, and the threshold considered as the success criterion. This threshold is defined by applying a relative error criterion of $5\%$, using the following formula:

$$\mathrm{target}=\mathrm{global}\mathrm{minimum}+0.05\times \mid \mathrm{global}\mathrm{minimum}\mid .$$

(22)

In the particular case of the Sum of Different Powers function, whose global minimum is zero, a target equal to ${10}^{-3}$ is arbitrarily set to avoid zero values.

Different families of parametric distributions (e.g., Gaussian, log-normal, etc.) are fitted to the parameters estimated using the Hilbert approximation implemented in the proposed method and compared with SBO. Ten runs of the code are performed, each with a randomly generated DoE. An optimization budget equal to 15 × dimension is set. The size of the DoE is set to 10 × problem dimension. The hyperparameter $\beta$ was assigned values between 10 and 100. Two different values are used for the parameter $\gamma$ :

• The first, $\gamma$, implemented within the optimization algorithm, takes values between 0.01 and 0.05.
• The second, ${\gamma }_{gb}=0.5$, is used to divide the DoE into two subsets (good and bad) for the construction of the apriori.

4.2. Benchmark Problems

In addition to the box plot presented in Figure 3(a)-3(f), we provide complementary numerical results in Tables 1-6. For each test function and each prior distribution used in the HSBO framework, we indicate:

• the maximum and minimum values reached by the function (corresponding to the 10 runs) obtained in the last iteration. This allows us to capture the worst and best performance of the algorithm.
• the mean of the function values at the last iteration (corresponding to the 10 runs) and the standard deviation of these values.
• the success rate, defined as the number of executions in which the function reached a value below the target.

These statistics provide additional information about the stability and robustness of our proposal and allow us to evaluate how often a method successfully identifies a near-optimal solution.

To evaluate the performance of the algorithm across all iterations in the ten runs, the simple regret metric was used. Following the definition in [29,30], simple regret at iteration t is expressed as:

$${\widehat{r}}_t^{+}=f\left({\mathbf{x}}^{*}\right)-\underset{{\mathbf{x}}_t\in {\mathcal{D}}_{1:t}}{min}f\left({\mathbf{x}}_t\right),$$

(23)

where ${\mathbf{x}}^{*}$ is the global optimum of the objective function and ${\mathcal{D}}_{1:t}$ denotes the set of points evaluated up to iteration t.

The mean simple regret obtained in the 10 runs is presented on a logarithmic scale in Figure A1 in Appendix B, which provides a clear comparison of the convergence speed of each configuration throughout the entire sequence of iterations.

4.2.1. Test Function a: Gramacy & Lee Function

Figure 3(a) shows the results obtained with the Gramacy & Lee function. It can be seen that both SBO and the proposed method with the Pearson distribution perform well. However, the mean of simple regret graph presented in Appendix B indicates that the Pearson distribution leads to faster convergence.

Table 1 shows that all methods manage to reach the theoretical minimum at least once. In terms of success rate, both SBO and the proposed method with Pearson distribution manage to find values below the target in all runs.

Table 1. Performance at Final Iteration ($f_{15}$) over 10 Runs: Gramacy & Lee 1D Function

Method	Min (${\mathbf{f}}_{15}$)	Max (${\mathbf{f}}_{15}$)	Mean (${\mathbf{f}}_{15}$)	Std(${\mathbf{f}}_{15}$)	Success rate ($\mathbf{f}\left(\mathbf{x}\right)<$ target)
SBO (baseline)	-0.8690	-0.8492	-0.8648	$7.6\times {10}^{-3}$	10/10
HSBO with normal prior	-0.8690	-0.6622	-0.8407	$6.5\times {10}^{-2}$	9/10
HSBO with log-normal prior	-0.8690	-0.6619	-0.8271	$8.7\times {10}^{-2}$	8/10
HSBO with exponential prior	-0.8690	-0.5185	-0.8099	$1.3\times {10}^{-1}$	8/10
HSBO with Rayleigh prior	-0.8690	-0.5185	-0.7922	$1.3\times {10}^{-1}$	7/10
HSBO with Pearson prior	-0.8690	-0.8492	-0.8648	$7.6\times {10}^{-3}$	10/10

Table 1. Performance at Final Iteration ($f_{15}$) over 10 Runs: Gramacy & Lee 1D Function

Method	Min (${\mathbf{f}}_{15}$)	Max (${\mathbf{f}}_{15}$)	Mean (${\mathbf{f}}_{15}$)	Std(${\mathbf{f}}_{15}$)	Success rate ($\mathbf{f}\left(\mathbf{x}\right)<$ target)
SBO (baseline)	-0.8690	-0.8492	-0.8648	$7.6\times {10}^{-3}$	10/10
HSBO with normal prior	-0.8690	-0.6622	-0.8407	$6.5\times {10}^{-2}$	9/10
HSBO with log-normal prior	-0.8690	-0.6619	-0.8271	$8.7\times {10}^{-2}$	8/10
HSBO with exponential prior	-0.8690	-0.5185	-0.8099	$1.3\times {10}^{-1}$	8/10
HSBO with Rayleigh prior	-0.8690	-0.5185	-0.7922	$1.3\times {10}^{-1}$	7/10
HSBO with Pearson prior	-0.8690	-0.8492	-0.8648	$7.6\times {10}^{-3}$	10/10

4.2.2. Test Function b: Cross-in-Tray 2D Function

Figure 3(b) shows the results obtained with the Cross-in-Tray function. It can be seen that both the normal and log-normal a priori distributions offer better performance compared to the SBO reference. Although SBO shares a low median, it has an outlier, indicating that in at least one run it deviated significantly from the optimum.

Table 2 confirms the observations made in the box plot: the normal and log-normal distributions achieve a success rate of 100% (10/10) and have a standard deviation lower than that of the baseline, indicating less variability between runs.

4.2.3. Test Function c: 2D Branin Function

Figure 3(c) shows the results obtained with the Branin function, where it can be seen that although the a priori-based configurations show competitive performance, in this particular case they fail to exceed the baseline (SBO).

Table 3 summarises the results obtained for the Branin function. It can be seen that the SBO method achieves the best average performance, with the lowest standard deviation and the highest success rate (6/10). In comparison, the HSBO configurations show greater variability and a reduced success rate, indicating that, in this case, the a priori assumptions do not provide a significant improvement over the baseline.

4.2.4. Test Function d: Hartmann 3D Function

Figure 3(d) shows the results for the Hartmann 3D function. It can be observed that both the baseline and the HSBO configurations with different prior distributions perform very similarly. The medians and dispersion of the relative errors are practically equivalent, suggesting that, in this case, none of the approaches offers a significant advantage over the others. This conclusion is reinforced by the numerical results presented in Table 4.

4.2.5. Test Function e: Sum of Different Powers (4D) Function

Figure 3(e) and Table 5 present the objective function values obtained at the end of the optimisation for the Sum of Different Powers function in dimension 4. Given that the global minimum is equal to zero, the most effective methods are those that achieve values closest to this threshold. In this regard, both SBO and HSBO with normal and Pearson distributions achieve outstanding performance, with very low final values and a 100% success rate. In contrast, configurations with lognormal, exponential and Rayleigh distributions have significantly higher function values and greater variability, indicating lower effectiveness for this problem.

4.2.6. Test Function f: Hartmann 6D Function

Finally, Figure 3(f) and Table 6 show the results obtained with the Hartmann 6D function. Although both SBO and the HSBO variants with different prior distributions achieve similar values close to each other, none of them comes close enough to the global minimum to meet the success criterion, which is reflected in a zero success rate for all cases. This result could indicate that, in problems of greater complexity or dimensionality, the number of iterations considered is not sufficient to allow effective exploration of the search space and reach the global optimum. Nevertheless, we highlight an interesting qualitative phenomenon in this particular case: the HSBO with lognormal prior continues to improve even after SBO has settled at a constant value. On average, SBO reaches a stable optimum around iteration 20. In contrast, HSBO-lognormal continues to improve around iteration 60, reaching a lower final average simple regret than SBO. This suggests an enhanced capability for exploitation for HSBO for this specific parametric prior.

4.3. Application: Shape Optimization of a Solid in Linear Elasticity Under Uniaxial Loading

A structural optimization problem is addressed in which the initial geometry, which is L-shaped, is parameterized using isogeometric analysis (IGA) on a finite element mesh constructed from NURBS surfaces [31]. The mechanical solver used is based on a MATLAB implementation of IGA in 2D linear elasticity, originally developed by Vinh Phu Nguyen (Johns Hopkins University) for educational and academic purposes1. This implementation includes the construction of stiffness matrices from NURBS, adaptive h refinement, and the treatment of natural and essential boundary conditions.

This code has been adapted to incorporate a shape optimization problem, in which the geometry of the L-shaped plate is perturbed by three design parameters that vertically modify certain NURBS control points. The design vector is denoted:

$$\mathbf{x}=\left[\begin{array}{c}{\epsilon }_{y1}\hfill \\ {\epsilon }_{y2}\hfill \\ {\epsilon }_{y3}\hfill \end{array}\right]\in \mathcal{X}\subset {\mathbb{R}}^3$$

where each variable ${\epsilon }_{yi}$ represents a vertical disturbance applied to a NURBS control point.

Figure 4. L-shape under uniaxial loading. Figure 4(a) shows the mechanical problem: the lower right edge of the structure is subjected to a uniform tensile load, while the left and lower edges are restrained in the x and y directions, respectively. Figure 4(b) shows how the control points that define the geometry are distributed. Figure 4(c) shows the initial configuration of the domain, and, Figure 4(d) illustrates the displacement in the x direction, represented by a contour plot.

Figure 4. L-shape under uniaxial loading. Figure 4(a) shows the mechanical problem: the lower right edge of the structure is subjected to a uniform tensile load, while the left and lower edges are restrained in the x and y directions, respectively. Figure 4(b) shows how the control points that define the geometry are distributed. Figure 4(c) shows the initial configuration of the domain, and, Figure 4(d) illustrates the displacement in the x direction, represented by a contour plot.

The goal is to minimize the total volume while ensuring that the average horizontal displacement matches a prescribed target. The optimization problem is therefore formulated as follows:

$$\underset{\mathbf{x}\in \mathcal{X}}{min}f\left(\mathbf{x}\right):=V\left(\mathbf{x}\right)+\lambda ·\left|{\overline{u}}_x\left(\mathbf{x}\right)-u_{\mathrm{target}}\right|$$

where:

• $V\left(\mathbf{x}\right)$ is the total volume of the generated structure,
• ${\overline{u}}_x\left(\mathbf{x}\right)={\displaystyle \frac{1}{|{\mathsf{\Gamma }}_D|}}{\sum }_{i\in {\mathsf{\Gamma }}_D}{u}_x^{\left(i\right)}$ is the average horizontal displacement on the right boundary ${\mathsf{\Gamma }}_D$,
• $u_{\mathrm{target}}=0,0013$ is the target displacement,
• $\lambda =100$ is a regularization coefficient.

Each evaluation of $f\left(\mathbf{x}\right)$ requires solving an isogeometric elasticity FE problem:

$$\mathbf{K}\left(\mathbf{x}\right)\mathbf{u}\left(\mathbf{x}\right)=\mathbf{F},$$

where $\mathbf{K}\left(\mathbf{x}\right)$ is the stiffness matrix depending on the geometry modified by $\mathbf{x}$, and $\mathbf{u}\left(\mathbf{x}\right)$ is the associated displacement field. This formulation guarantees mechanical consistency via the imposed equilibrium, while the volume term, although approximate, plays a regularizing role, guiding the solution towards more compact forms without explicitly restricting the space of admissible designs.

To initiate the optimization process, an initial DoE comprising 12 samples was generated using Latin hypercube sampling (LHS). The admissible set $\mathcal{X}$ is defined by the lower and upper bounds ${[-0.50.5]}^3$. To evaluate the acquisition function at each iteration, a Monte Carlo sampling of ${10}^6$ points was extracted from the design space. The hyperparameters of the proposed optimization algorithm were established as follows: the quantile threshold $\gamma =0.01$ was used in the definition of the probability of improvement, and the scale coefficient $\beta =100$. Furthermore, to divide the DoE into ‘good’ and ‘bad’ subsets, ${\gamma }_{gb}=0.5$ response values were taken to construct a priori distributions representing the promising and unpromising regions of the design space.

Figure 5 shows the evolution of the volume throughout the optimization process. The initial volume corresponds to the original L-shaped design. The first 20 iterations were performed using the proposed HSBO method. During this phase, the volume was progressively reduced, especially in the first iterations. Once these 20 iterations were completed, the best solution obtained with HSBO, as well as the trained Gaussian surrogate, were taken as the starting point and surrogate objective functions for a local gradient-based optimization procedure, implemented using Matlab’s fmincon function. The algorithm used by fmincon (by default, interior-point) was configured with standard tolerance settings. This second stage successfully converged towards the optimum, achieving a total volume reduction of approximately 30% compared to the initial design. In terms of computational cost, HSBO used 12 functional evaluations to build the DoE, followed by 20 more evaluations during the iterative process, for a total of 32 evaluations. Since fmincon was applied to the Gaussian Process surrogate, no additional functional evaluation was required during this refinement step. Overall, the entire process involved 32 functional evaluations. This hybrid approach illustrates HSBO’s versatility and its efficiency in training the Gaussian surrogate with scarce data and with high accuracy around the region surrounding the optimum.

Figure 6 shows the results of the optimized L. Figure 6(a) shows the final positions of the three control points that act as design variables. Their vertical displacements generate the diagonal notch responsible for the volume reduction. Figure 6(b) shows the resulting refined NURBS mesh, which has uniform quality elements and a deformation field without local concentrations, confirming the mechanical consistency of the optimized shape. Figure 6(c) shows the contour of the horizontal displacement $u_x$, which exhibits a nearly linear gradient between the embedded edge and the right edge; the uniform color band on ${\mathsf{\Gamma }}_D$ confirms that the average displacement reaches the target $u_{\mathrm{target}}=1.3\times {10}^{-3}\mathrm{m}$ with minimal dispersion.

5. Discussion

In this article, we have presented HSBO, a version of BO that combines Bayesian-Hilbertian UQ modeling techniques within the framework of Bayesian optimisation. The role of HA is the construction of structured performance priors over the search space, identifying regions of promising performance to guide the Bayesian enrichment process. We have performed several synthetic tests to evaluate the behavior of this method, revealing competitive performance in selected cases.

The numerical tests, carried out on various test functions with varying characteristics and dimensions, reveal that the proposed method provides a systematic framework for prior construction that can enhance classical BO performance when problem structure aligns with the selected parametric family. The improved performance can manifest as increased success rate, faster convergence, or closer proximity to the optimum on a fixed budget, with respect to SBO. Even in situations where the global minimum is not identified, the method converges towards regions very close to the optimum, which represents an excellent starting point for gradient-based refinement methods.

This latter property has been particularly useful in the structural shape optimization application presented, where the aim was to minimize the volume of a solid, initially an L-shaped plate. The Bayesian-Hilbertian approach made it possible to identify, with only 20 iterations, a restricted region of the search space enclosing the optimal solution. This region, and the associated metamodel trained with the selected points, served as the initial basis for a gradient-based algorithm to refine the search, yielding an accurate solution with no additional evaluation of the underlying expensive objective function.

This approach is particularly relevant for engineering problems where expert priors are difficult to define, yet initial design data can be obtained with reasonable effort. In such cases, constructing a small exploratory design of experiments (DoE) is worthwhile, as it enables a principled, data-driven initialization of the optimization process with minimal reliance on human bias.

Future lines of inquiry include: the investigation of the impact of explicit modeling of dependencies between data when constructing the priors; sensitivity analysis of the hyperparametes involved, particularly $\gamma$ and $\beta$; consideration of other surrogates. Additionally, two promising complements of this work would include the manipulation of the priors to integrate constrained optimization into this framework, and guided biasing of the priors to boost convergence speed.