Hierarchical Gaussian Process-Based Bayesian Optimization for Materials Discovery in High Entropy Alloy Spaces

In materials discovery applications, GP and MTGP models offer powerful frameworks for predictive modeling and optimization. GPs are a robust tool for modeling objective functions in multi-objective BO, providing a probabilistic approach to regression that quantifies uncertainty in predictions [66]. Traditional BO algorithms often model each objective function independently using separate GPs, thus assuming no correlation between them. This independent modeling approach can overlook potential interrelationships among objective functions. In contrast, the MTGP approach models all objective functions jointly, capturing any correlations and leveraging shared information to improve predictive accuracy [67].

Consider a GP constructed for an objective function $f(x)$. The prediction at an unobserved location $\mathbf{x}$ given $N$ previously observed data points denoted by $\{\mathbf{X}_{N},\mathbf{y}_{N}\}$, where $\mathbf{X}_{N}=(x_{1},x_{2},\ldots,x_{N})$ and $\mathbf{y}_{N}=(f(x_{1}),f(x_{2}),\ldots,f(x_{N}))$, is given as:

where the mean $\mu(x)$ and variance $\sigma^{2}(x)$ are defined as:

Here, $k$ is the kernel function, $\mathbf{K}(\mathbf{X}_{N},\mathbf{X}_{N})$ is an $N\times N$ matrix where the $(m,n)$-th entry is $k(x_{m},x_{n})$, and $\mathbf{K}(\mathbf{X}_{N},x)$ is an $N\times 1$ vector with the $m$-th entry being $k(x_{m},x)$. The term $\sigma^{2}$ represents experimental noise. The kernel function captures the correlation between observations based on their relative distances [68]—these distances may not necessarily be defined in a Euclidean space. A common choice for the kernel function is the squared exponential kernel:

where $d$ is the dimensionality of the input space, and $l_{h}$ is the characteristic length-scale determining the strength of correlation in the $h$-th dimension.

The MTGP model extends this by jointly modeling multiple objective functions, taking into account their inter-task correlations [69]. To achieve this, we define a task-specific covariance matrix $\mathbf{K}_{f}$ to capture the similarities between tasks. The overall covariance function for the MTGP model is constructed using the Kronecker product of the task covariance matrix $\mathbf{K}_{f}$ and the input covariance matrix $\mathbf{K}_{x}$:

Here, $\otimes$ denotes the Kronecker product. The variable $\mathbf{y}_{N}$ in the equations is a column vector, where each task’s data is ordered in columns and stacked vertically.

This MTGP model allows us to capture the correlations between different objective functions, thereby improving the accuracy and efficiency of the optimization framework. By leveraging the shared information between tasks, the MTGP can enhance the performance of the design optimization process, leading to more accurate predictions and faster convergence to optimal solutions [70]. Inference in the MTGP model follows the standard GP equations for the mean and variance of the predictive distribution, extended to the multi-task setting. The mean prediction for a new data point $x^{*}$ for task $l$ is given by:

where $\mathbf{k}_{f}^{l}$ selects the $l$-th column of $\mathbf{K}_{f}$, $\mathbf{k}_{x}^{*}$ is the vector of covariances between the test point $x^{*}$ and the training points, and $\Sigma$ is the combined covariance matrix given by $\Sigma=\mathbf{K}_{f}\otimes\mathbf{K}_{x}+\mathbf{D}\otimes\mathbf{I}$, where $\mathbf{D}$ is a diagonal matrix representing noise variances. Here, $y$ represents the observed data across all tasks, stacked appropriately. The use of MTGP allows for the transfer of information between tasks, leveraging the inter-task correlations to improve predictive performance, especially when data is sparse for some tasks. This joint modeling approach is particularly beneficial in materials discovery, where the properties of materials are often correlated, and exploiting these correlations can lead to more efficient and accurate predictions [71, 72].

The hyperparameters of the MTGP model, including the parameters of the kernel functions and the task covariance matrix, are learned by maximizing the marginal likelihood of the observed data [66]. This can be done using gradient-based optimization methods, ensuring the positive semi-definiteness of $\mathbf{K}_{f}$ through parametrization such as the Cholesky decomposition [69]. Overall, the MTGP model provides a flexible and powerful approach for multi-objective optimization, capturing inter-task correlations and enhancing the efficiency of the design framework.

Deep Gaussian Processes (DGPs) extend the concept of GP to hierarchical models, providing a more powerful and flexible approach to modeling complex data [73]. A DGP consists of multiple layers of latent variables, where each layer is modeled by a GP. This hierarchical structure allows for the representation of more abstract features at higher layers, which can improve predictive performance in tasks with complex underlying structures [74].

In a standard GP, we model a set of training input-output pairs, $\{\mathbf{X},\mathbf{Y}\}$, using a latent function $f(x)$ drawn from a GP prior, and we infer the distribution of $f$ given the data. Observed data points are denoted by $\{\mathbf{X},\mathbf{Y}\}$, where $\mathbf{X}=(x_{1},x_{2},\ldots,x_{N})$ and $\mathbf{Y}=(y_{1},y_{2},\ldots,y_{N})$. In the case of the DGP, the inputs to the GP at each layer are governed by the outputs of the GP from the previous layer. Formally, if we consider a DGP with $H$ hidden layers, the generative process can be described as:

where $\mathbf{y}_{n}$ are the observed data points, $\mathbf{h}_{n}^{i}$ are the latent variables at the $i$-th layer, and $\mathbf{x}_{n}$ are the discrete data points used as the inputs for the highest layer in the hierarchy. Each function $f^{i}$ is modeled as a GP with an appropriate covariance function [75].

Training a DGP involves marginalizing over the latent variables in each layer, which is analytically intractable. To address this, we use a variational approach to approximate the marginal likelihood [76]. The variational distribution $Q$ is introduced to approximate the true posterior distribution of the latent variables. The variational lower bound on the marginal likelihood is given by:

where $Q$ is factorized as:

and $\mathbf{U}^{i}$ are inducing variables for the $i$-th layer, which is a smaller set of data points designed to approximate the posterior for each layer while managing computational complexity. The terms in the variational lower bound can be computed analytically, allowing for efficient optimization of the model parameters and variational distributions [77].

This Bayesian training approach not only provides a principled way to learn the model parameters but also enables automatic determination of the appropriate structure for the deep hierarchy through the use of automatic relevance determination (ARD) priors [78]. The ARD priors help in identifying the most relevant dimensions at each layer, facilitating effective model complexity control.

Overall, DGPs offer a robust framework for capturing complex dependencies in data, making them particularly suitable for applications in materials discovery where the relationships between variables can be highly nonlinear and intricate [73, 79].

In BO, acquisition functions guide the search for the optimal solution by balancing exploration and exploitation. Two popular acquisition functions are the Upper Confidence Bound (UCB) and Expected Hypervolume Improvement (EHVI).

In BO with DGPs, the distinction between isotopic and heterotopic sample sites(query sites or data points) is crucial for understanding the operational differences between Isotopic DGP-BO (iDGP-BO) and Heterotopic DGP-BO (hDGP-BO)[69].

Isotopic data refers to scenarios where sample sites are shared across all tasks. This means that the same input locations are used for every task, simplifying the modeling process since the correlation structure between tasks can be more easily captured due to the shared sampling sites. This is the approach taken by iDGP-BO in this work, which assumes that the tasks are closely related and can be effectively modeled using shared sample points or query sites. While this assumption reduces complexity and improves computational efficiency, it may miss task-specific nuances that arise from differing sample sites, potentially leading to suboptimal performance in complex, real-world scenarios where such assumptions do not hold. Mathematically it can be described as,

where $x_{i}$ represents the shared sample locations, and $f_{1}(x_{i})$ and $f_{2}(x_{i})$ represent the outputs for the respective tasks at those locations.

Heterotopic data, on the other hand, allows query sites to differ between tasks, providing greater flexibility in modeling. This is the foundation of hDGP-BO , where each task has its own set of sampling sites. This approach leverages the additional flexibility provided by different query locations to better capture task-specific correlations, making it possible to exploit these for more effective optimization. By allowing for different sample sites, hDGP-BO can perform more nuanced and accurate predictions, particularly in cases where the tasks exhibit significant variability in their data distributions. This added complexity requires more sophisticated modeling but can result in superior performance, as evidenced in our experiments. It can be expressed as,

where $x_{i,1}$ and $x_{i,2}$ represent different sample locations for each task.

iDGP-BO uses a DGP that assumes all tasks share the same sampling sites. Here, all the queries across all the task are always made on same query locations based on EHVI acquisition function. This isotopic approach simplifies the correlation structure between tasks, making it easier to model but potentially less flexible in handling task-specific variations. The advantage of iDGP-BO lies in its simplicity and computational efficiency, which can be beneficial in scenarios where tasks are expected to be closely related.

hDGP-BO, however, uses a more complex framework that allows each task to have its own sampling sites (heterotopic approach). This model leverages DGP to capture more intricate correlations between tasks, which can be particularly beneficial when tasks are less related or exhibit significant individual characteristics. hDGP-BO also incorporates an acquisition strategy where less expensive tasks can be queried more frequently, using techniques like UCB to reduce uncertainty, thereby informing the sampling strategy for more expensive tasks. This makes hDGP-BO not only more flexible but also potentially more accurate in scenarios where task-specific data is critical. In our strategy, hDGP-BO makes queries using same input points one-third of the times using EHVI. During other steps, it queries only BM and volume at different input locations, based on the highest UCB of respective tasks. In this way hDGP-BO can explore or be informed about CTE at a certain location by only querying the BM or volume at that specific location. This occurs because of the correlated structure of DGP. Thus, knowing one task at a specific input location informs the GP about associated tasks at that location without explicitly querying those associated tasks.

In the workflow of this current work, hDGP-BO queries CTE, BM and volume 500 times using EHVI during the whole BO procedure (isotopic queries). In addition to that, hDGP-BO queries only BM and volume twice using UCB before each isotopic queries. So, in total hDGP-BO acquires CTE data 500 times and BM, volume data 1500 times. Whereas, models executing only isotopic queries (cGP-BO, MTGP-BO & iDGP-BO) acquires CTE, BM and volume data 500 times each. hDGP-BO leverages this additional relatively inexpensive BM and volume data to carry out the optimization task in a more efficient and robust manner.

As the underlying simulation platform for our BO, we employed the pyiron workflow framework [88] to conduct high-throughput screening simulations for sampling the FeCrNiCoCu high entropy alloy (HEA) space. pyiron provides workflows to calculated material properties independent of the simulation engine, including thermal expansion coefficeints (CTE), equilibrium bulk modulus (BM) and equilibrium volume (V). We calculate those with an atomic composition resolution of 3%, amounting to a total of roughly 60,000 configuration to sample the composition space. For benchmarking the BO schemes the target properties were calculated for all compositions a priori. Subsequently, the BO schemes queried from this consolidated database for optimization, which accelerates benchmarking the different BO schemes and facilitates systematic performance analysis.

We emphasize that, in a real experimental design scenario, the properties would only be evaluated at the query points in an on-the-fly fashion, a mode of operation that our workflow is designed to support. The pyiron workflow framework facilitates the seamless integration and automation of various simulation tools and workflows, significantly enhancing the efficiency and scalability of high-throughput screenings. Conceptually, pyiron functions as a high-level task-based interface to the simulation engine responsible for the material property calculations. In our case, the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) [89, 90] is employed to efficiently evaluate material properties across a large number of configurations. Additionally, pyiron provides a user-friendly Python interface for setting up, executing, storing, and analyzing high-throughput molecular dynamics (MD) simulations at large computational scales on high-performance computing (HPC) resources, which is critical for the coupling with Bayesian Optimization (BO).

To efficiently sample the composition space and ensure a realistic representation of the random nature of HEAs, we utilized Special Quasirandom Structures (SQS) [91, 92] with 32 atoms, generated with the sqsgenerator package  [93]. SQS are designed to mimic the perfectly random atomic arrangements of alloys in finite supercells, providing more accurate simulation results. The atomic configuration were evaluted using the embedded atom method (EAM) interatomic potential developed by Deluigi et al., which was designed for simulating FeNiCrCoCu [94]. That being said, as the purpose of this study is to quantify the performance of various multi-task BO approaches, the absolute accuracy of the physical predictions obtained from the interatomic potential is not critical, as the general statistical properties of the correlations within and between tasks is expected to control performance.

For our workflow, each simulation began with rescaling the simulation cell based on the concentration-averaged nearest neighbor distance of its end members. This is followed by a structure optimization including relaxation of the simulation cell volume and internal degrees of freedom, resulting in an equilibrated structure at 0K. For calculation of equilibrium BM, these equilibrated structures are strained from -15% to +15% with an additional equilibration of the internal degrees of freedom. The BM is then calculated as the second derivative of these energy volume curves. Starting from the same equilibrated structure at 0K the CTE is determined by varying the simulations temperature from 15K to 1500K in 5K intervals over a 6ps molecular dynamics trajectory in the NPT ensemble [95, 96] with a 1fs time step. By recording the change of volume with increasing temperature, the CTE is calculated providing insights into the thermal behavior of the alloys [97]. These MD simulations are carried out in $3\times 3\times 3$ supercells of the 32 atoms SQS cells, resulting in 864 atoms cells to ensure statistical accuracy [98, 99].The calculation of CTE differs from that of BM. The CTE require calculations of volume(along with energy minimization and NPT ensemble trajectories) at 297 different temperatures. Whereas, for calculation of equilibrium BM, there’s no dependence on temperature unlike CTE. Thus, the cost associated with BM calculations is lower than that of CTE from a computational perspective. The CTE calculations requires 29800 more force evaluations (within a NPT ensemble), compared to the cost of calculation in case of BM.

In total 3.5 billion atomic configurations were evaluated to construct the consolidated database for the roughly 60,000 compositions with over 750,000 individual simulation tasks being recorded in the pyiron database. This highlights the need of a workflow framework like pyiron for systematically evaluating the HEA composition space.

The main objective of this work is to compare the performances of cGP-BO, MTGP-BO, iDGP-BO and hDGP-BO in two specific multi-objective optimization tasks, concerning FeCrNiCoCu HEA design space. The first optimization task is to minimize CTE and maximize BM. The second optimization task involves maximization of both CTE and BM in the vast HEA space. Both of the objectives has their respective engineering applications, which will be discussed in details afterwards.

Firstly, the results of our high-throughput atomistics simulations for the FeCrNiCoCu alloy system are presented using Uniform Manifold Approximation and Projection (UMAP) representations, as shown in Figure 2 [100, 101]. UMAP, a dimension reduction technique, was chosen for its ability to preserve both local and global structure in high-dimensional data, making it particularly suited for visualizing the complex relationships in our alloy composition space.

UMAP offers several advantages over traditional dimension reduction methods [101]:

• 1.

  It preserves topological structure across multiple scales, allowing for the visualization of both fine-grained local relationships and broader global patterns.

• 2.

  It can handle non-linear relationships in the data, which are common in materials science datasets.

• 3.

  It is computationally efficient, enabling the visualization of large datasets typical in high-throughput studies.

In our analysis, UMAP transforms the five-dimensional compositional space of the FeCrNiCoCu alloy system into a two-dimensional representation. Each point in these visualizations corresponds to a unique alloy composition, with color gradients indicating the magnitude of different material properties. This approach provides an intuitive understanding of how properties vary across the compositional landscape, allowing for rapid identification of trends, clusters, and potential high-performance regions. These plots effectively capture the complex interplay between composition and properties, offering insights that might be obscured in traditional analysis methods. The use of UMAP allows us to visually explore correlations between different properties and identify compositional regions of interest for further investigation.

Figure 2 reports the properties of the FeCrNiCoCu alloy space in UMAP scheme, calculated using the workflow described in the last section. Figure 2(a) reports the mapping of the compositional space of the FeCrNiCoCu alloy system onto a UMAP. Colors indicate region where the alloy contains an element at more than 50% concentration, while no single element reaches this threshold in grey regions. Pure materials are found at the vertices of the projection. Figure 2(b) illustrates the distribution of thermal expansion coefficients across different alloy compositions in this same projection, showing significant variation and potential candidates with low CTE. Figure 2(c) shows the bulk moduli for the same set of alloys, which is crucial for identifying compositions with high mechanical strength. Finally, Figure 2(d) depicts the equilibrium volumes, providing additional context for the structural properties of the alloys. Note that this projection does not guarantee continuity of the various properties in the 2D plane. It should therefore be taken primarily as a visualization aid.

In addition to having UMAP plots for the properties individually, attempts have been made to understand the inter-relationship between the properties. The pairwise plots in Figure 3 and the correlation table in Table 1 provide a deeper understanding of the inter-dependencies between thermal expansion coefficients, bulk moduli, and volumes. Table 1 reports a statistical correlation analysis, including Pearson correlation coefficients, Spearman coefficients, and mutual information between CTE, BM, and volume. The Pearson and Spearman coefficients provide insights into the linear and monotonic relationships between the properties, respectively, while mutual information captures both linear and non-linear dependencies [102].

The weak correlations observed between these properties suggest the possibility of complex non-linear relationships [102, 103]. It is important to note that while weak linear correlations can indicate non-linear relationships, they are not definitive proof [104]. Nonetheless, the presence of weak correlations in our high-dimensional materials data motivates the exploration of more sophisticated analysis techniques [105]. Thus, the potential for using multi-task learning models to exploit these potentially complex non-linear relationships is a motivating factor for our subsequent analysis [67, 106, 107].

These correlations arise from the underlying physical mechanisms governing the behavior of HEAs, such as the complex interplay between atomic size mismatch, lattice distortion, and chemistry  [108, 109]. Empirical models to predict the thermal expansion based on the equilibrium volume, equilibrium bulk modulus and its derivative like proposed by Moruzzi et al. [110] fail for these HEAs, which highlights the need for improved regression surrogates capable of learning complex relationship like MTGP or DGP for HEAs. Finally, leveraging these correlations can lead to more efficient exploration and optimization of the HEA design space.

In this section, we characterize and compare the relative efficiency of four multi-task BO approaches: cGP-BO, MTGP-BO, iDGP-BO and hDGP-BO. We apply these approaches to multi-objective material design tasks within the FeCrNiCoCu HEA space, focusing on the simultaneous minimization of CTE and maximization of BM (also, maximization of both as a case study in the next subsection), while using equilibrium volume as an auxiliary task.

This optimization scenario is particularly relevant for advanced technological applications requiring materials with exceptional dimensional stability under varying thermal conditions while maintaining high mechanical strength. Such materials are crucial in industries like aerospace, where components must withstand extreme temperature fluctuations without compromising structural integrity. For instance, in the design of next-generation aircraft engines or hypersonic vehicles, alloys with low CTE and high BM are essential for maintaining precise geometries and resisting deformation under high-stress, high-temperature environments [108, 111].

Our investigation aims to demonstrate how advanced BO techniques, particularly MTGP-BO and DGP-BO, can leverage inter-task correlations and complex, non-linear relationships in the material property space to accelerate the discovery of HEAs with these desired characteristics. By comparing these methods against cGP-BO, we seek to quantify the potential improvements in optimization efficiency and highlight the advantages of multi-task learning approaches in navigating the vast and complex design space of HEAs. Among all the different variants of BO, our novel hDGP-BO has demonstrated the most significant improvements over all other BO variants. The primary objective in this scenario was to minimize the CTE while maximizing the BM.

Figure 4 reports the evolution of the hypervolume during an iterative BO procedure based on different models, where an iteration is defined by a single query to all tasks using EHVI, except for the heterotopic DGP-BO (hDGP-BO) model where one iteration of hDGP-BO comprises of two steps of queries to BM and volume using UCB, and one step of query to CTE, BM and volume using EHVI. So, during 500 iteration, hDGP-BO acquired 500 CTE, 1500 BM and 1500 volume data points, whereas, cGP-BO, MTGP-BO & iDGP-BO acquired 500 datapoints of CTE, BM and volume each. One of the objectives is to investigate whether hDGP-BO can use these extra to accelerate the mapping of the Pareto front with respect to isotopic BO variants.

Figure 5 & Figure 6 visualizes the queries in the objective space and design space for the first 100 queries, and the Pareto optimality metric across the query space, respectively. The Pareto optimality metric is defined by the euclidean distance to the nearest Pareto point, scaled between 0 and 1, where 1 implies being farthest from the nearest Pareto optimal points and 0 indicate the Pareto set itself. The KDE (kernel density estimation) contounrs represent the distribution of queries made by 20 independent replications of the whole BO loop.

Figure 4 along with Figure 5 & Figure 6 compare the performance and querying behavior of cGP-BO, MTGP-BO, isotopic DGP-BO (iDGP-BO), and heterotopic DGP-BO (hDGP-BO). DGP-BO and MTGP-BO models, in general, exhibit a more thorough exploration of the objective space, efficiently identifying regions with low CTE and high BM. The superior performance of MTGP-BO and DGP-BO can be attributed to their ability to exploit correlations between different objectives, as noted by prior works [112, 70, 113]. This leads to fewer iterations to identify optimal solutions, thereby reducing computational costs and time.

A careful examination of Figure 4 along with Figure 5 & Figure 6 further highlights differences between iDGP-BO, hDGP-BO & MTGP-BO. The hDGP-BO was by far performing best, showing excellent hypervolume improvements compared to others, covering more than 95% of the actual maximum hypervolume within the first 50 iterations. In comparison, MTGP-BO (the second best model) required more than 100 iterations to achieve the same coverage. The performance bounds obtained in multiple replications of hDGP-BO was also narrower, demonstrating a lower sensitivity to initial conditions, presumably due to the additional information gained from the additional interleaved queries of BM and volume using UCB.

This interpretation is supported by the fact that iDGP-BO failed to match the performance of hDGP-BO, while sharing the same architecture. iDGP-BO was also performing worse than MTGP-BO in latter stages. This suggests that UCB-based uncertainty-minimization improves the performance of DGP. The query pattern in Figure 4 & Figure 5 indicates that iDGP-BO made queries roughly randomly, showing no clear sign of systematic exploration. MTGP-BO was able to perform better than iDGP-BO and ultimately reached the hypervolume attained by hDGP-BO at latter stages. But, it was not able to achieve that as fast as hDGP-BO. Figure 4 & Figure 5 indicates the querying pattern of MTGP-BO is less concentrated around the Pareto front compared to hDGP-BO, leading to lower performance.

This phenomenon may be related to the simpler linear relationship assumed between tasks by the MTGP’s structure.

This more efficient exploration and exploitation of the design space by hDGP-BO corroborates findings by Hernandez-Lobato et al. (2016) [114], who demonstrated the effectiveness of advanced BO models in multi-objective optimization tasks through the use of acquisition functions based on information-theoretic approaches. Our obsevations indicate that the combination of an UCB heterotopic sampling strategy and of the non-linear cross-task representation lead to the superior performance of hDGP-BO.

In the second optimization scenario, the objective was to maximize both the CTE and the BM. Finding an alloy with both high BM and high CTE is generally challenging due to the typical inverse relationship between these properties in most materials [115], since i) still interatomic interactions tend to oppose thermal expansion [116, 117] and an empirically-observed inverse correlation between anharmonicity (which is the driving force for CTE) and bond strength [118, 119].

Materials possessing both high BM and high CTE simultaneously, while relatively rare, find valuable applications in specific engineering contexts. In thermal management systems for electronics, they can provide structural stability while allowing for thermal expansion to match surrounding components. This property combination is particularly useful in electronic packaging, where managing thermal stresses is crucial [120]. High-temperature sensors operating in extreme temperature environments also benefit from materials that maintain structural integrity while being sensitive to temperature changes. The high bulk modulus ensures the sensor’s structural stability under pressure, while the high CTE allows for enhanced sensitivity to temperature variations [121]. In the field of thermal energy storage, materials that can withstand high pressures and exhibit significant volume changes with temperature are desirable. These properties are particularly relevant for phase change materials and other advanced thermal storage systems, where the ability to maintain structural integrity under varying thermal conditions is critical [122]. Micro-electro-mechanical systems (MEMS) and certain types of actuators also benefit from materials combining high stiffness with high thermal responsiveness. The high bulk modulus provides the necessary structural rigidity, while the high CTE allows for thermally-driven actuation. This combination is particularly valuable in the development of advanced MEMS devices and smart materials for various applications, including aerospace and robotics [123]. It’s important to note that the specific requirements for ”high” values of these properties can vary depending on the application context. The development of materials that optimize both properties simultaneously remains an active area of research in materials science and engineering, driving innovation in various technological fields.

Following the structure of the previous section, Figure 7 reports the evolution of the hypervolume during an iterative BO procedure based on the four different models. Figure 8 & Figure 9 visualizes the queries in the objective space and design space for the first 100 queries. Figure 8 also displays the Pareto optimal points for reference, whereas Figure 9 reports the corresponding Pareto optimally metric, defined as described above.

Figure 7 along with Figure 8 & Figure 9 displays the difference in performance and querying behavior of cGP-BO, MTGP-BO, iDGP-BO, and hDGP-BO in this specific scenario. As in the previous problem, DGP-BO and MTGP-BO outperform cGP-BO, generating a more systemic exploration of the objective space, efficiently identifying regions with high CTE and high BM. Figure 7 indicates that MTGP-BO and both variants of DGP-BO (heterotopic and isotropic) outperform cGP-BO in terms of hypervolume improvement per iteration.

A closer scrutiny of Figure 7 along with Figure 8 & Figure 9 reveals that, unlike the first case, iDGP-BO no longer performs worse than MTGP-BO, but instead performs better, on par with hDGP-BO. The difference in performance between both DGP methods on the one hand and MTGP-BO and cGP-BO on the other hand is now even wider compared to the previous scenario (c.f., Figure 7 vs Figure 4 ) Both DGP methods systematically hone in to the optimal region, while the queries from MTGP-BO & cGP-BO are more diffuse and less systematic, leading to a poor exploration of Pareto-optimal solutions, which is consistent with the poor performance exhibited in Figure 7.

The superior exploration and search capability of hDGP-BO can be attributed to the same reasons as stated in the previous scenario. But, improvement of iDGP-BO performance over MTGP-BO along with overall decrease in performance of cGP-BO and MTGP-BO points to the difficulty of the current optimization goals. It can be inferred that, as the problem became harder, the simple linear assumption of MTGP-BO falters during earlier iterations. It demonstrated slower convergence to maximum attainable hypervolume compared to iDGP-BO and hDGP-BO. But, at latter stages, when MTGP-BO was able to explore enough data points, it slowly achieved the same level of coverage as that of DGP based BOs. So, at latter stages MTGP-BO matches the performance with that of DGP based BOs. But, during earlier iterations (within the first 50 iterations) DGP based BOs showed exceptional performance in covering close to 90% of the maximal hypervolume, compared to MTGP-BO which achieved only 75% coverage after 100 iterations. Additionally, cGP-BO performs even worse due to the enhanced difficulty of this problem. The difficulty faced by MTGP-BO and cGP-BO is also reflected in the querying pattern( Figure 8 & Figure 9). In both of the cases, the hDGP-BO proposed in this paper is the most robust.

To further analyze the differences in performance among multiple variants of BO, an analysis based on mutual information(MI) is carried out. The prediction of the four BOs were recorded in each iteration at 5000 points(sampled through dirchlet distribution) for each tasks. Then, the mutual information between the two objectives were calculated from the predicted 5000 points. This was done at each iteration.

In the first scenario, Figure 10(a) shows that, only the predicted MI of hDGP-BO gets close to the true MI. Whereas, iDGP-BO struggles to achieve any significant MI due to its unsystematic query pattern(Figure 6(c)). MTGP-BO display high predicted MI due to the fact that they assume linear correlation between tasks. High MI translates into simpler and linear relationship. The cGP-BO also shows high MI without having correlated kernel.It’s because the different tasks start out with same prior, making the tasks highly correlated at unexplored regions. The agreement between predicted and true MI translates into the performance difference(Figure 4). The BO having the best agreement(hDGP-BO) shows the best performance.

The second scenario(Figure 10(b)) reaffirms the hypothesis ”ability to achieve correct predicted MI translates into better optimization performance”. Here, it is observed that iDGP-BO is able to capture the correct predicted MI(unlike first scenario). As a result, the iDGP-BO also performs better(Figure 7), unlike the first case. Moreover, the predicted MI from cGP-BO and MTGP-BO shows significant deviation from true MI in this second case. As a result, their optimization performance also deteriorates compared to the first scenario(Figure 7).

In this work, predicted MI per iteration has been demonstrated to be useful in quantifying performance for optimization problems having multiple correlated tasks. Further research is warranted to investigate the use of predicted MI in different multi-objective optimization problems having multiple correlated tasks. It can be proposed as a metric of the degree to which a model can exploit correlations.

In this subsection, we analyze the impact of adding random Gaussian field noise to the Bayesian optimization process. The addition of noise is essential for evaluating the robustness and stability of the optimization algorithms under real-world conditions, where measurements and simulations are often subject to various sources of uncertainty [124, 125].

Gaussian field noise is characterized by its mean and covariance structure, which can be defined as:

where $k_{\epsilon}(\mathbf{x},\mathbf{x}^{\prime})$ is the covariance function of the noise, typically chosen to be a squared exponential kernel:

Here, $\sigma_{\epsilon}^{2}$ represents the noise variance, and $l_{\epsilon}$ is the length scale of the noise [125]. In this work, The lengthscale is chosen to be 0.20. The variance is selected such that, the third standard deviation does not exceed more than half of the spread for each tasks. As a result, CTE and BM have variances of 20 and 100 respectively. This implies uncertainty in CTE data having 4.5 $\times$$10^{-6}$ K^-1 standard deviation and in BM having 10 GPa standard deviation.

To assess the influence of this noise on the optimization performance, we conducted experiments by adding random Gaussian field noise to the objective functions during the Bayesian optimization process. The performance of the optimization algorithms was evaluated both before and after the addition of noise [126].

Figure 11 and Figure 12 illustrate the hypervolume improvement per iteration for the two optimization scenarios: minimization of TEC with maximization of BM, and maximization of both TEC and BM, respectively. The figures compare the performance of the optimization algorithms before and after the addition of random Gaussian field noise.

The results indicate that the addition of noise generally leads to a decrease in the optimization performance advantages between the models. The DGP-BO models still show better performance compared to the cGP-BO & MTGP-BO model. But, the performance advantages of DGP-BOs decrease. The inherent uncertainty between multiple replications also increases. This happens due to the fact that, addition of gaussian field noise obscures the complex non-linear relationship between tasks. Thus, it becomes harder for DGPs or MTGPs to take advantage of the inter-relationship between tasks. This further proves that the DGPs and MTGPs performs better then cGP-BO by exploiting correlations. Otherwise, obscuring the correlation through noise addition would not have affected the performance difference between them.

The hDGP-BO still demonstrates robustness to obscured correlation between tasks. This robustness can be attributed to the sophisticated kernel structure and uncertainty reduction steps using UCB. In the context of HEA design, the presence of noise is a common challenge due to the complexity of the systems and the limitations of experimental and computational methods [111]. The robustness of hDGP-BO to noise makes them particularly suitable for handling the uncertainties inherent in HEA property data, ensuring more reliable optimization results. The incorporation of random Gaussian field noise in the optimization process serves as a valuable test for evaluating the robustness of Bayesian optimization algorithms. The findings underscore the advantages of using advanced models hDGP-BO in scenarios with inherent uncertainties, which are common in real-world materials discovery tasks [127].