Dynamic Exploration-Exploitation Trade-Off in Active Learning Regression with Bayesian Hierarchical Modeling

In machine learning, exploration refers to generating novel information in the search space while exploitation focuses on improving decisions (e.g., maximizing reward) based on the available information. Although this is a generally accepted definition of exploration and exploitation [44], its usage varies depending on the context. Particularly, in the regression setting, exploration refers to sampling data from the unobserved regions in the search space to gather more (potentially novel) information about the black-box function, such as local minima/maxima. In contrast, exploitation aims at accurately capturing the black-box function in the regions where it is highly unpredictable and has a sharp change or discontinuity. As mentioned in the foregoing, it is possible to encode exploration and exploitation separately in an acquisition function [34], however, it is challenging to dynamically balance between the two [45]. We consider a general framework for sampling data points that maximizes a linear combination of acquisition functions aimed at simultaneous exploration and exploitation as presented in Equation (1). The next data is queried through pure exploitation when $\eta=0$, and pure exploration when $\eta=1$. By increasing $\eta$ from $0$ to $1$, the degree of exploration increases in the constructed acquisition function.

To estimate $\eta$ dynamically, we consider two levels of variability. First is the variability in queried data emerging from the noise and measurement errors in $y$ at each querying stage. The second level of variability is associated with the dynamic nature of $\eta$ between querying stages due to the nature of the unknown black-box function (e.g., jumps or discontinuities). This bi-level nature of variability motivates the need for a hierarchical model. A Bayesian hierarchical model allows us to capture this hierarchy while simultaneously encoding the uncertainty at each of the levels. The trade-off parameter $\eta_{j}$ captures the first level of variability at querying stage $j$ given the set of hyper-parameters $\bm{\theta}$. We model $\{\eta_{j}|\bm{\theta}_{j}\}\stackrel{{\scriptstyle iid}}{{\sim}}p(\eta|{\bm{\theta}_{j}})$. Here, the conditional independence of $\{\eta_{j}|\bm{\theta}_{j}\}$ tells us that the trade-off parameters at one sampling stage are exchangeable (de Finetti’s theorem [46]), but not independent. To capture the variability in $\eta$ across multiple sampling stages, we model $\{\bm{\theta}_{j}|\bm{\psi}\}\stackrel{{\scriptstyle iid}}{{\sim}}p(\bm{\theta}|\bm{\psi})$ where $\bm{\psi}$ is another set of hyper-parameters fixed from a priori assumptions. De Finetti’s theorem from this conditional independence shows that the trade-off parameters at subsequent querying stages are exchangeable as well. Indeed, this aligns with our assumption that the trade-off parameter at any querying stage is not completely independent of the previous actions (exploration or exploitation). In the next two subsections, we discuss the prior and the sampling distributions employed in this study.

In the absence of any known functional form of $f(\mathbf{x})$, we consider Gaussian process regression (GPR) as our underlying learning model. GPR considers a distribution over the underlying function $f(\mathbf{x})$ and aims to specify the black-box function by its mean and covariance. Due to the noise inherited in labeled data without the learner’s knowledge, the output $y$ comprises $f(\mathbf{x})$ as $y=f(\mathbf{x})+\varepsilon$ where $\varepsilon\sim\mathcal{N}\left(\mathbf{0},\sigma_{n}^{2}\right)$. Observing the noisy outputs, ${\mathbf{y}}_{o}$, GPR attempts to reconstruct the underlying function $f(\mathbf{x})$ by removing the contaminating noise, $\varepsilon$ [53]. We denote the queried dataset as $({X}_{o},\mathbf{y}_{o})$. GPR imposes a zero-mean Gaussian process prior over the noisy outputs such that,

where $X_{u}$ is the set of unlabeled testing points. The posterior distribution at the test data is given as $\{\hat{\mathbf{f}}|X_{o},\mathbf{y}_{o},X_{u}\}\sim\mathcal{N}\left(\bar{\mathbf{f}},\text{cov}(\hat{\mathbf{f}})\right)$ where,

In this section, we present some of the most commonly used acquisition functions and characterize them either as an exploration or an exploitation strategy based on the literature and our extensive simulation studies.

$\bullet\quad\textbf{Improved Greedy Sampling: }$ As mentioned in Section 2, the concept of improved greedy sampling (iGS) was introduced to ensure diversity in the queried data [10]. At every iteration, iGS determines the unexplored region by searching over the entire region in both the input and output spaces. It queries the next data which is located the farthest from its nearest training point or in an unobserved region according to the then prediction by the regression model. The acquisition function at $\mathbf{x}_{j}$ is calculated as $\min\left(u(\mathbf{x}_{j})v(\mathbf{x}_{j})\right)$ where $u(\mathbf{x}_{j})=\|\mathbf{x}_{j}-X_{o}\|_{2}$ and $v(\mathbf{x}_{j})=\|\bar{f}_{j}-\mathbf{y}_{o}\|_{2}$ refer to the distance of $\mathbf{x}_{j}$ with training points in input and output space respectively. Here $\mathbf{y}_{o}$ is the set of the observed output or the labels at training points $X_{o}$, $\bar{f}_{j}$ is the predicted output at test point $\mathbf{x}_{j}$, and $\|\cdot\|_{2}$ is the L2 norm. IGS defines the distance at $\mathbf{x}_{j}$ as $u(\mathbf{x}_{j})v(\mathbf{x}_{j})$ instead of $u(\mathbf{x}_{j})+v(\mathbf{x}_{j})$ or $\left(u(\mathbf{x}_{j})\right)^{2}+\left(v(\mathbf{x}_{j})\right)^{2}$ due to the possibility of significantly different scale of $u$ and $v$ which can hamper the latter two formulas with the dominance of one measure over another. It then selects the next data, $\mathbf{x}_{\text{iGS}}^{*}$, such that,

$\bullet\quad\textbf{Query by Committee: }$ Successfully applied in classification and regression-based problems, Query by Committee (QBC) is a framework rooted in the concept of utilizing an ensemble of hypotheses [31]. Maintaining a committee of models, QBC queries the data where the committee members disagree the most about a measure of criteria. The committee members are all trained on the same set of training points but with competing hypotheses or regression models. Considering a committee of $Q$ models denoted as $h_{1},h_{2},\ldots,h_{Q}$, we define the measure of disagreement between two models $h_{l}$ and $h_{p}$ at $\mathbf{x}_{j}$ as the absolute difference of prediction by the respective models at that point, or $|h_{l}(\mathbf{x}_{j})-h_{p}(\mathbf{x}_{j})|$. Then the acquisition function at $\mathbf{x}_{j}$ is defined as $\max_{l,p}(|h_{l}(\mathbf{x}_{j})-h_{p}(\mathbf{x}_{j})|)$, representing the maximum disagreement at $\mathbf{x}_{j}$ where $l,p=1,2,\ldots,Q$. Then the next data, $\mathbf{x}_{\text{QBC}}^{*}$ is selected such that,

$\bullet\quad\textbf{Maximum variance: }$ This strategy defines the acquisition function as the variance predicted by the regression model[54]. A learner’s expected error can be decomposed into noise, bias, and variance [8]. Since the noise is independent of the model or data, and bias is invariant given a fixed model, minimizing the variance is intuitively guaranteed to minimize the future generalization error of the model [5]. Equation (10) provides the Gaussian process posterior covariance, and the diagonal elements of $\text{cov}(\hat{\mathbf{f}})$ provide the predicted variance, $\mathbb{V}$. The next data is queried following,

$\bullet\quad\textbf{Maximum entropy: }$ Another active learning strategy formulated based on uncertainty is the maximum entropy strategy [5]. Entropy is an information-theoretic measure that often has been defined as the amount of information needed to “encode” a distribution and has been related to the uncertainty of the underlying model [5]. Thus minimizing model entropy can reasonably lead to revealing the model uncertainty. Shannon’s entropy has been used as an acquisition function named maximum entropy sampling [7] or uncertainty sampling [5]. Since the posterior prediction of the Gaussian process follows a multivariate normal distribution, Shannon’s entropy of the distribution,

As touched on briefly in the introduction, the existing approaches have tried achieving the trade-off between exploration and exploitation using different methodologies [38, 55, 56, 16]. Here we have listed a few that we will compare BHEEM with.

$\bullet\quad\textbf{Static trade-off: }$ Many of the previous studies have held on to static trade-offs between exploration and exploitation [38]. Prior information can be used to decide on the trade-off (e.g., equal importance to both or more importance to one based on the nature of the function). Existing approaches have conducted trial and error with different trade-off values between exploration and exploitation between appropriate exploration and exploitation acquisition functions to conduct the simulated experiments, and then select the one that promises overall better accuracy [55]. By fixing the exploration-exploitation trade-off throughout the whole learning process, this method reduces the computation complexity to a great extent. Nevertheless, the same trade-off cannot be expected to demonstrate similar performance for every function (as can be seen in the experimental results).

$\bullet\quad\textbf{Probabilistic trade-off: }$ Inspired by simulated annealing [57] and built on the $\epsilon-$decreasing greedy algorithm [56], this algorithm updates the exploration-exploitation combination in a probabilistic approach [16]. For instance, the exploration probability is defined as $p_{R}=\alpha^{t-1}$ where $\alpha$ is less than 1 and $t$ is the current time step or iteration number. To query new data, a uniform random variable $Z$ is generated, if $Z\leq p_{R}$, we consider $\eta=1$, and exploration is performed, otherwise, we consider $\eta=0$, and exploitation is applied [16]. Hence the strategy starts with pure exploration. The exploration probability decays gradually, and the learning model gets more prone to pure exploitation over time. Intuitively, this transition is practical since the initial queried data should focus more on exploration, and after learning the overall trend of the function, we can identify the irregular regions via exploitation. However, the transition rate depends on the choice of $\alpha$ which we fixed at $0.7$ following Elreedy et al. [16]. But again, the transition rate should depend on the nature of the function and fixing it will decrease the efficiency of the learning process.

To demonstrate the efficacy of BHEEM, we employ the Matérn 3/2 kernel function in the GPR model to avoid too rough (exponential kernel) or too wavy (squared exponential kernel) prediction. The kernel function is defined as,

While applying the QBC strategy, we maintained a committee of ten Gaussian process models, each with a different kernel function; (i) squared exponential, (ii) exponential, (iii) Matérn 3/2, (iv) Matérn 5/2, (v) rational quadratic, (vi) product of dot product and constant kernel, (vii) product of i and iii, (viii) product of i and iv, (ix) product of ii and iii, (x) product of ii and iv [53]. For the generation of the data, we fixed the signal-to-noise (SN) ratio to 10 which we define as the decibel of the ratio of signal power to the power of the data noise. For the performance evaluation, we calculated the root mean squared error as,

To demonstrate the performance of BHEEM against the existing strategies, we considered the six following functions.

• •

  $F_{1}(x)=\begin{cases}3.5\exp\left(-\frac{(x-10)^{2}}{200}\right)+\epsilon,&\qquad\quad\qquad\quad\quad\text{if $x\leq 25$}\\ 8-3.5\exp\left(-\frac{(x-35)^{2}}{200}\right)+\epsilon,&\qquad\qquad\quad\quad\quad\text{otherwise}.\end{cases}\qquad\qquad\;\;\quad\qquad\quad(17)$

• •

  $F_{2}(x)=\sin(x)+2\exp(-30x^{2})+\epsilon\qquad\qquad\qquad\quad\quad\quad\;\;x\in[-2,2]\qquad\qquad\quad\qquad\;\;\quad(18)$

• •

  Three-hump camel function:
  $F_{3}(\mathbf{x})=2x_{1}^{2}-1.05x_{1}^{4}+x_{1}^{6}/6+x_{1}x_{2}+x_{2}^{2}+\epsilon\qquad\qquad\quad x_{1},x_{2}\in[-5,5]\qquad\qquad\qquad\;\;(19)$

• •

  Six-hump camel function:
  $F_{4}(\mathbf{x})=(4-2.1x_{1}^{2}+x_{1}^{4}/3)x_{1}^{2}+x_{1}x_{2}+(4x_{2}^{4}-4)x_{2}^{2}+\epsilon\quad x_{1},x_{2}\in[-2,2]\qquad\qquad\qquad\;(20)$

• •

  Hartmann 3D function:
  $F_{5}(\mathbf{x})=-\sum_{i=1}^{4}\alpha_{i}\exp{\left(-\sum_{j=1}^{3}A_{ij}(x_{j}-P_{ij})^{2}\right)}+\epsilon\quad\quad x_{i}\in[0,1]\;\forall i=1,2,3\quad\qquad\quad(21)$
  

  where $\alpha=(1.0,1.2,3.0,3.2)^{T}$, $A=\begin{bmatrix}3.0&10&30\\ 0.1&10&35\\ 3.0&10&30\\ 0.1&10&35\end{bmatrix}$, and $P=10^{-4}\begin{bmatrix}3689&1170&2673\\ 4699&4387&7470\\ 1091&8732&5547\\ 381&5743&8828\end{bmatrix}$

• •

  $F_{6}(\mathbf{x})=\frac{1}{2}\sum_{i=1}^{10}x_{i}^{2}-\cos\left(2\pi x_{i}\right)+\epsilon\quad\qquad\qquad\qquad\quad\quad x_{i}\in[-5,5]\;\forall i=1,2,\ldots,10\quad\;(22)$

First, we implement BHEEM by employing iGS and QBC as the pure exploration and exploitation strategy respectively. With three random initial queried data, Figure 3 shows the result of iGS, QBC, and BHEEM after the $12^{\text{th}}$ iteration for the underlying function from Equation (18). Unsurprisingly, iGS spreads out its queried data in both the $x$ and $y$ direction predicting the smooth portions of the function competently, but fails to fit the peak of the function due to the sudden change near $x=0$. QBC exploits this uncertainty in predictions across different committee members who lead to different predictions near the sharp peak. Therefore, QBC queries most of the data in this region while ignoring other regions in the process. BHEEM queries data near $x=0$ enough times to fit the peak there but has also explored and provided satisfactory prediction in other regions achieving an overall lower error than the other two acquisition functions.

Continuing to employ iGS and QBC as pure exploration and exploitation strategies respectively, Figure 4 (a-f) compares the root mean square error (RMSE) of different methodologies while learning Equation (17-22) and summarises the result using boxplots. For each of the examples, we have set the upper limit on the total number of experiments to be 100. The key observation in the boxplots is that our proposed methodology, BHEEM, consistently achieves a lower generalization error with a smaller number of queried data than any other employed strategy. Between iGS and QBC, we observe the latter to be more efficient in some cases (Figure 4(c)), and iGS in others (Figure 4(e)).

Among the one-dimensional problems, BHEEM achieved the lowest RMSE compared to others after the query of the first $10$ data in the case of $F_{1}(x)$ as observed in Figure 4(a). IGS was the closest one in this case and achieved only about $5\%$ higher RMSE on an average. While predicting $F_{2}(x)$, BHEEM converged the fastest and performed the best as shown in Figure 4 (b). The second best strategy was iGS which scored about $60\%$ higher RMSE on average from BHEEM after querying $25$ data. Among the two-dimensional problems, in $F_{3}(\mathbf{x})$, BHEEM surpassed others most of the time and at the $50^{\text{th}}$ addition of point, achieved about $8.5\%$ lower RMSE from QBC which scored the closest to BHEEM for the corresponding functions. For $F_{4}(\mathbf{x})$, the consistency of lower generalization error and fast convergence persisted for BHEEM. For the three-dimensional $F_{5}(\mathbf{x})$, BHEEM and QBC were the two strategies achieving the lowest RMSE across the querying stages. BHEEM was able to achieve significantly lower RMSE than the pure exploration strategy, iGS. In the ten-dimensional problem of $F_{6}(\mathbf{x})$ (Figure 4(f)), all the active learning strategies are observed to compete with each other over the learning process. The pure exploration strategy, iGS, is performing better than all other acquisition functions at the $50^{\text{th}}$ addition of data whereas QBC seems to achieve the lowest RMSE at the $100^{\text{th}}$ addition. In most of these cases, all active learning strategies performed better than the random sampling strategy. Our result was affected by a few decisions including our choice of kernel, sampling algorithm and proposal distribution, number of points in the initial set of queried data, etc. However, the overall result demonstrates that our proposed approach tends to perform at least tantamount to pure exploration or pure exploitation.

To compare the computational time taken by the applied methodologies, Figure 4 also presents the average time for each methodology to query 100 data during the learning process. Our proposed methodology, BHEEM, takes a significantly larger processing time, about $2.5$ and $2$ times higher than the duration of iGS and QBC, respectively. Hence, BHEEM is more suitable for offline applications that involve time-consuming physical experiments than the onlineones, which is the case for most applications in manufacturing and materials that involve conducting time consuming physical experiments.

Figure 5 presents the percentage improvement of BHEEM from the pure exploration (iGS) and pure exploitation (QBC) calculated from the average RMSE obtained in the simulated examples. Corresponding to each function, the four bars of each color represent the percentage improvement from exploration using iGS (blue) and exploitation using QBC (yellow) after querying the fifth, tenth, $25^{\text{th}}$, and $50^{\text{th}}$ data. Here, the positive and negative bars indicate increased and decreased accuracy of the BHEEM respectively, compared to pure exploration and pure exploitation. We observe that among the twenty-four cases of displayed results, only one showed deterioration where BHEEM had lower accuracy than both the pure strategies. In all other stages, BHEEM had significant improvement and achieved lower RMSE compared to at least one of them while surpassing result of both pure strategies in sixteen of them.

We also compare our proposed methodology with static and probabilistic updates of $\eta$. For the static trade-off, we considered $\eta=0.25,0.5,0.75$ during the learning process individually and calculated the RMSE for each of them. To note, $\eta=0$ and $\eta=1$ refer to pure exploitation and pure exploration respectively which we have already compared with BHEEM in Figure 5 and Figure 6. For the probabilistic update of $\eta$, we followed the strategy described in Section 3.3.2. A heatmap for the average RMSE scaled for each of the comparative trade-off methodologies is presented in Figure 6 where blue and red cells represent more and less accurate models respectively. The heatmap clearly shows that no one value of $\eta$ works well for every function, or even for every iteration in one function. BHEEM distinctly provides better results than all the static and probabilistic trade-offs.

Finally, Figure 7 presents the average improvement achieved by BHEEM from pure exploration and exploitation strategy across the functions. Here, we have used iGS and maximum variance as the exploration strategy in Figure 7(a) and Figure 7(b) respectively, and QBC as the exploitation strategy in both cases. In Figure 7(a), we observe about $7\%$ and $21\%$ lower average RMSE from pure exploration (iGS) and exploitation (QBC) respectively. In Figure 7(b), we observe about $5.7\%$ and $2.3\%$ lower average RMSE from pure exploration (maximum variance) and exploitation (QBC) respectively. Overall, our proposed methodology of trade-off always promises better accuracy than pure exploration and exploitation irrespective of the choice of strategies.

The layered ternary carbides and nitrides with the general formula $M_{n+1}AX_{n}$ are called MAX phase materials where $M$ and $A$ are early transition metal and A-group elements respectively, whereas $X$ is Nitrogen or Carbon and $n$ is an integer between 1 and 4 [58]. Their layered structures kink and delaminate the materials during deformation resulting in an unusual and unique combination of both ceramic and metallic properties which makes them attractive candidates for structural and fuel coating applications. In this study, we intended to predict the lattice constants of MAX phase materials by analyzing their compositions and elastic constants using the data from Aryal et al. [59]. The lattice constant is an important piece of information to define the overall lattice structure, which helps model the microstructure evolution. To represent the discrete categorical composition features into numeric input to the models, we used one-hot encoding, a popular approach replacing the categorical variable with as many variables as categories [60]. Assuming that we wish to differentiate between two materials with the same elements in M and A, and the same value of $n$, one of them is a carbide, while the other is a nitride. It is possible to represent this information with two binary variables using one-hot encoding where $0$ and $1$ represent the non-existence and existence of that element in the material respectively.

Each category value is represented as a 2-dimensional, sparse vector of $1$ for one of the dimensions, and $0$ for the other. In general, for variables of cardinality $d$, the one-hot encoding would transfer it to $d$ number of binary variables where each observation indicates the presence (1) or absence (0) of the dichotomous binary variable [60]. In the MAX phase problem, we have ten elements in M, twelve elements in A, and two elements in X. Hence, after using one-hot encoding and adding the numerical variable for n, we have a total of 25 variables representing the composition of the materials. Including the composition and the elastic constants, there is a total of 30 predictors to predict the lattice constant. Figure 8 provides a comparative analysis of RMSE of different methodologies employed in this case study. Our proposed methodology (combining iGS and QBC) outperformed the other strategies and surpassed at least one of the pure exploration or exploitation in almost every iteration. BHEEM achieved about $10.4\%$ and $11.3\%$ average improvement from the pure exploration and exploitation strategy respectively.

In this section, we discuss the convergence of sampled $\eta$ at one querying stage, its relation with the choice of hyperparameters, and the dynamics $\eta$ during the active learning process. The results in this section are all based on the learning process of $F_{1}(x)$.

$\bullet\quad\textbf{Convergence of Metropolis: }$ We visualize the progression of the Metropolis algorithm in Figure 9 for a representative simulation example of $F_{1}(x)$. It plots the accepted $\hat{\eta}$ at one stage of the querying process. The chains were initiated with different ranges of values for $\hat{\eta}$ and continued till 10,000 samples were selected. The green, orange, and blue lines represent MCMC chains initiated with sampled values of $\eta$ ranging between $(0,0.2)$, $(0.4,0.6)$, and $(0.8,1)$ respectively. For each range of initial values, we plotted only three MCMC chains for better visualization. Figure 9(a) and Figure 9(b) demonstrate the chains till 50 and 10,000 accepted samples respectively. It is evident from the figures that irrespective of the initial sample, MCMC chains at one stage of learning converge to a close neighborhood of $\eta$ very quickly.

$\bullet\quad\textbf{Dynamics of $\mathbf{\eta}$: }$ Progression of the posterior mean of $\eta$ during successive query stages while learning $F_{1}(x)$ is presented in Figure 9(c). Dotted line is used only for visualization of the sample path and has no physical meaning. The values of $\eta$ closer to $1$ and $0$ encourages exploration and exploitation respectively. Every other value of $\eta$ combines exploration and exploitation in a specific ratio. Based on this representative sample path, we first note that BHEEM is never stuck with exploration or exploitation which is the case when the trade-off parameter is pre-defined. A close investigation of the sampled $\eta$ shows that BHEEM continues to either explore, exploit, or combine both towards optimally learning the function as more data points are queried. Note that the sample path of the mean $\eta$ depends on the noise level therefore may be different on each run.

$\bullet\quad\textbf{Sensitivity of hyperparameters: }$ Through the sensitivity analysis of the standard deviation of proposal distribution, $\tau$, and threshold parameter, $\nu$, visualized in Figure 10, we will discuss the dependence of convergence of RMSE and selection of $\eta$ on both the hyperparameters. We repeated our experiments by changing $\nu$ as $0.0001,0.001,0.01,0.1,$ and $0.25$ while keeping $\tau$ fixed. Boxplot of root mean square error (RMSE) across different querying stages in Figure 10(a) indicates that different values of $\nu$ do not have a significant effect on RMSE the since we execute the MCMC sampling long enough to ensure convergence. However, the noise in data compels us to query the data differently each time, and the value of $\eta$ at the previous stages affect $\eta$ in subsequent stages. Therefore, selected values of $\eta$, presented as a boxplot in Figure 10(b), have a large variance across the learning process. Due to the near indifference to the result, we have chosen $\nu=0.001$ for all experimentation. Figure 10(c) and Figure 10(d) present the RMSE and the $\eta$ selected respectively for different values of $\tau$, $0.001,0.01,0.1,$ and $0.25$, while keeping $\nu=0.001$. Here also, we observe the behavior of RMSE and $\eta$ to be indifferent towards the value of $\tau$. But due to our previous discussion on the relation between $\tau$ and the speed of convergence of Metropolis algorithm in the last paragraph of Section 3.1.2, we chose $\tau=0.1$ throughout all experiments.