Exploratory Landscape Validation for Bayesian Optimization Algorithms

Abstract

Bayesian optimization algorithms are widely used for solving problems with a high computational complexity in terms of objective function evaluation. The efficiency of Bayesian optimization is strongly dependent on the quality of the surrogate models of an objective function, which are built and refined at each iteration. The quality of surrogate models, and hence the performance of an optimization algorithm, can be greatly improved by selecting the appropriate hyperparameter values of the approximation algorithm. The common approach to finding good hyperparameter values for each iteration of Bayesian optimization is to build surrogate models with different hyperparameter values and choose the best one based on some estimation of the approximation error, for example, a cross-validation score. Building multiple surrogate models for each iteration of Bayesian optimization is computationally demanding and significantly increases the time required to solve an optimization problem. This paper suggests a new approach, called exploratory landscape validation, to find good hyperparameter values with less computational effort. Exploratory landscape validation metrics can be used to predict the best hyperparameter values, which can improve both the quality of the solutions found by Bayesian optimization and the time needed to solve problems.

1. Introduction

Continuous optimization problems with high computational complexity in terms of objective functions arise in many fields of engineering, material design, automatic machine learning, and others [1]. In real-world optimization problems, evaluating the value of an objective function requires time-consuming computational experiments or simulations. The time needed to solve an optimization problem is primarily defined by the total number of objective function evaluations performed by the optimization algorithm. When solving computationally expensive optimization problems, the main termination criterion for the algorithm is the maximum number of objective evaluations, which is usually called the computational budget [2].

Problems with a limited computational budget are best solved by optimization algorithms that utilize surrogate models of the objective function [2]. The most promising candidates for objective evaluation are identified by refining and exploring surrogate models at every optimization iteration. Hence, both the quality of surrogate models and the parameters of the model exploration procedure can have an impact on the optimization algorithm’s efficiency. Many studies have focused on tuning the parameters of the model exploration procedure, specifically the type and the parameters of the acquisition function, as well as the parameters of the algorithm for optimizing it [3,4]. This article focuses on ways to improve the efficiency of optimization algorithms by increasing the quality of surrogate models.

Surrogate models are usually built by approximation algorithms with tunable hyperparameters. Numerous studies have been conducted to find the best hyperparameter values that maximize the accuracy of surrogate models and, hence, the optimization algorithm’s efficiency [5,6,7]. The process of finding such hyperparameter values is called hyperparameter tuning. Here are some examples of hyperparameters being tuned for different types of surrogate models:

• Degree of polynomial regression models;
• Size of the hidden layer for models based on neural networks;
• Kernel type and parameters for models based on Gaussian process (GP) regression;
• Number of trees and tree depth for models based on random forest regression.

This article considers the tuning of kernel parameters for models based on GP regression, which includes, for example, the scale mixture parameter $\alpha$ of the rational quadratic kernel or the $\nu$ parameter of the Matern kernel [8].

In order to perform hyperparameter tuning, the hyperparameter efficiency metric is defined based on the quality estimate of a surrogate model built with a given hyperparameter vector. The model quality is usually measured on a test sample of points using some error approximation metric, such as mean absolute error, mean squared error, $R^2$ score, etc. [9,10]. In cases of a strict computational budget, maintaining a separate test sample by evaluating a costly objective function is not sufficient. To increase the reliability of hyperparameter efficiency estimates, resampling procedures are used, which involve the building of multiple models for each hyperparameter vector, for example, the well-known cross-validation procedure [10]. Using cross-validation for hyperparameter efficiency estimation significantly increases the computational cost of hyperparameter tuning. This article introduces a new approach for estimating hyperparameter efficiency called exploratory landscape validation, along with new efficiency metrics that require less computational effort, and that, in some cases, outperform cross-validation in terms of solution quality. One of the metrics is based on a so-called ranking preservation indicator [11], calculated on an extended training sample. Another metric is evaluated by comparing the variability maps of an objective function [12] constructed for a training sample and a surrogate model.

Much of the recent research on the parameter tuning of optimization algorithms, which also includes the hyperparameter tuning of surrogate models, has focused on using exploratory landscape analysis (ELA) algorithms to estimate the characteristic features of optimization problems [13,14,15]. The hyperparameter prediction process is based on the following idea. Given the ELA feature vectors, which incorporate the specifics of problems, and the best hyperparameter values found for those problems, a machine learning (ML) algorithm is used to build a tuning model to predict the best hyperparameter values or the best approximation algorithm [12,16,17]. According to such an ELA-ML approach, the tuning model is then used to identify hyperparameter values that are reasonably effective when solving similar problems based on the feature vectors of those problems. In [16], for example, the authors constructed tuning models that are based on different classification algorithms to predict the best surrogate modelling algorithm using the ELA features of the problem. Despite the fact that additional computational effort is required to collect the training data and build a tuning model, the increase in efficiency of the optimization algorithm outweighs the computational costs in a long-term perspective. This article also explores the benefits of using tuning models that are built using the proposed metrics for finding effective hyperparameter values.

The rest of the article is organized in the following way. Section 2 presents a statement for a global continuous optimization problem, describes the canonical form of the Bayesian optimization algorithm, and formalizes the hyperparameter tuning and prediction problems. In Section 3, a new approach to the quality estimation of surrogate models, called landscape validation, is introduced, and includes a variety of criteria for selecting and predicting the best hyperparameter values. Section 4 starts with the general setup for computational experiments, followed by a description of the experiments performed with the hyperparameter tuning and hyperparameter prediction approaches based on the proposed hyperparameter efficiency metrics; then, the experimental results are discussed.

2. Bayesian Optimization with Hyperparameter Tuning and Prediction

The following statement of the continuous global optimization problem $q$ is considered:

$$\underset{X\in D_X\subset {ℝ}^{\left|X\right|}}{min}f\left(X\right)=f\left(X^{*}\right)=f^{*},$$

(1)

where $X=\left(x_1,\mathrm{...},x_{\left|X\right|}\right)$ is the vector of continuous variables of size $\left|X\right|$; $f\left(X\right)\in {ℝ}^1$is a scalar objective function; $D_X$ is a convex region of the permissible variables’ values; $X^{*}$ is the vector of the variables’ values for which the objective function has the minimal value $f^{*}$. It is assumed that the region $D_X$ is formed by the set of inequalities $x_i^{-}\le x_i\le x_i^{+}, i\in \left(1:\left|X\right|\right)$, where $x_i^{-}$ and $x_i^{+}$ are the lower and upper bounds of the $i$-th variable, respectively. The problem (1) is referred to as the base optimization problem.

Due to the high computational complexity of the objective function, the total number of objective evaluations allowed for solving the problem (1) is limited to $N_{max}$, which is called the computational budget of the problem. When solving computationally expensive base problems, common practice is to build and explore surrogate models $\widehat{f}\left(X\right)$ of the objective function $f\left(X\right)$. By using a surrogate model $\widehat{f}\left(X\right)$, promising areas of the search region $D_X$ can be located approximately without utilizing the budget $N_{max}$. The algorithm for building surrogate models usually has tunable hyperparameters, the vector of the values of which is denoted as $P=\left(p_1,\mathrm{...},p_{\left|P\right|}\right)$.

In this section, the Bayesian optimization algorithm is described and the problem statements for hyperparameter tuning and hyperparameter prediction for the surrogate modeling algorithm are formulated. The aim of this section is to outline and formalize common approaches for solving computationally expensive optimization problems, on the basis of which, new approaches to the hyperparameter efficiency estimation will be presented in the next section.

2.1. Bayesian Optimization Algorithm

In contrast to other surrogate-based optimization algorithms, Bayesian optimization algorithms use surrogate models that are based on Gaussian process (GP) approximation. In GP models, each point $X$ of the search region $D_X$ is associated with a normal distribution of the predicted objective values $\mathcal{N}\left({\mu }_{\widehat{f}}\left(X\right),{\sigma }_{\widehat{f}}^2\left(X\right)\right)$, where ${\mu }_{\widehat{f}}$ is the mean and ${\sigma }_{\widehat{f}}$ is the standard deviation of $\widehat{f}$ values. Having the probability distribution over $\widehat{f}$ allows for selection of promising points for the objective function evaluation based on both predicted objective values and uncertainty estimates of those predictions. Formally, promising point selection is defined as the maximization problem of a so-called acquisition function. Hence, the tunable parameters of a Bayesian optimization algorithm include the type and parameters of the acquisition function, as well as the hyperparameters of the GP approximation algorithm, which are the type and parameters of a covariance function, also referred to as a kernel [5].

The Bayesian optimization algorithm includes the following steps (Figure 1):

• Generate the initial sample $L^0=\left\{\left(X^i,f^i\right), i\in \left(1:N_{init}\right)\right\}$ for training the first GP model, where $N_{init}<N_{max}$ is the initial sample size. Points $X^i\in D_X$ are chosen either randomly or according to one of the designs of the experiment algorithm, e.g., the Latin hypercube sampling (LHS) algorithm [18];
• Perform iterations $r\in \left(1:N_{iter}\right)$, where $N_{iter}=N_{max}-N_{init}$ as follows:
  • Build the surrogate model ${\widehat{f}}^r\left(X\right)$ using the current training sample $L^r$ and the vector $P$ of the hyperparameter values;
  • Select the next point $X^{i+1}$ for the objective function evaluation by optimizing an acquisition function, e.g., by minimizing the lower confidence bound (LCB) [19] function as follows:

    $$X^{i+1}=\underset{X\in D_X\subset {ℝ}^{\left|X\right|}}{\arg \min }{\mu }_{\widehat{f}}\left(X\right)-\kappa {\sigma }_{\widehat{f}}^2\left(X\right),$$

    (2)

    where $\kappa$ is a tunable parameter;
  • Evaluate the objective function for point $X^{i+1}$ to obtain the corresponding value $f^{i+1}$ and extend the training sample $L^{r+1}=L^r{\displaystyle \cup }\left(X^{i+1},f^{i+1}\right)$.
• The point that has the minimal corresponding objective value $\left(X^{*},f^{*}\right)\in L^{N_{iter}}$ is considered as the solution of the base problem (1).

Although the formal optimality of the best-found point $X^{*}$ is not guaranteed in any sense, the same notation as in the problem definition (1) is used for simplicity. The rest of the article focuses on ways to improve the optimization algorithm’s efficiency by tuning the GP hyperparameters.

2.2. Hyperparameter Tuning for a Bayesian Optimization Algorithm

The efficiency of Bayesian optimization can be improved by selecting, at each iteration, the vector $P^{*}$ that is the best for the training sample $L^r$ according to the hyperparameter efficiency metric $\varphi \left(L^r,P\right)$. Given the set of allowed hyperparameter values $D_P\subset {ℝ}^{\left|P\right|}$, the best vector $P^{*}$ is found by solving the following hyperparameter optimization problem at step 2a (Figure 2):

$$\underset{P\in D_P\subset {ℝ}^{\left|P\right|}}{opt}\varphi \left(L^r,P\right)=\varphi \left(L^r,P^{*}\right),$$

(3)

The hyperparameter efficiency metric $\varphi \left(L^r,P\right)$ is commonly defined based on approximation accuracy metrics, e.g., mean squared error, on a test sample.

In most cases, the size of sample $L$ in Bayesian optimization is already not sufficient for the search space dimension $\left|X\right|$ due to a very limited budget $N_{max}$. Since it would not be practical either to split sample $L$ into train and test parts or to spend the budget $N_{max}$ on collecting and updating a separate test sample, the cross-validation procedure is used to estimate the average accuracy of the surrogate models built with the given vector $P$. The corresponding hyperparameter efficiency metric ${\varphi }_{CV}\left(L^r,P\right)$ is calculated as follows:

$${\varphi }_{CV}\left(L^r,P\right)=\frac{1}{K}{\displaystyle \sum }_{k=1}^K{R}^2\left(L_k^r,P\right)=\frac{1}{K}{\displaystyle \sum }_{k=1}^K\left[1-\frac{{\displaystyle \sum }{\left(f_k^i-{\widehat{f}}_k^i\right)}^2}{{\displaystyle \sum }{\left(f_k^i-{\overline{f}}_k\right)}^2}\right],$$

(4)

where $K$ is the total number of cross-validation folds, $R^2\left(L_k^r,P\right)$ is the coefficient of determination, $L_k^r\subset L^r$ is the test sample for the $k$-th fold, $f_k^i$ and ${\widehat{f}}_k^i$ are the known and predicted objective values at the $k$-th fold correspondingly, and ${\overline{f}}_k$ is the mean of the known objective values at the $k$-th fold.

Using the metric ${\varphi }_{CV}\left(L^r,P\right)$ when solving the problem (3) requires building as many GP models for each vector $P$ as there are cross-validation iterations. Since building a GP model has $O\left({\left|L\right|}^3\right)$ complexity, where $\left|L\right|$ is the training sample size, the hyperparameter tuning process may become time consuming and unprofitable. In this article, it is proposed to develop new hyperparameter efficiency metrics that reduce the computational complexity of solving the problem (3) while maintaining the accuracy of solutions comparable to the metric ${\varphi }_{CV}\left(L^r,P\right)$.

2.3. Hyperparameter Prediction for a Bayesian Optimization Algorithm

Solving the problem (3) from scratch at each iteration of the Bayesian optimization is a straightforward but time-consuming approach to hyperparameter tuning. Modern approaches to parameter tuning or the selection of optimization algorithms involve combining ELA and ML to predict the most efficient algorithm or parameters for solving the base problem (1) [16,17]. The ELA-ML approach relies on the assumption that the best optimization algorithm or parameter values for solving similar optimization problems will also be nearly identical. The similarity of optimization problems can be estimated by a similarity measure between the corresponding ELA feature vectors.

Although many of the known ELA-ML approaches are developed for the best optimization algorithm selection, the same idea can be adapted for hyperparameter tuning the following way. The process of solving the problem (3) is divided into separate phases as is shown in Figure 3.

Exploration phase—before solving the base problem (1):

• Define a representative set $Q=\left\{q_i,i\in \left(1:M_Q\right)\right\}$ of test optimization problems $q_i$, where $M_Q$ is the number of test problems;
• Generate random training samples $\left\{L_{i,j}, j\in \left[1:M_L\right]\right\}$ of different sizes $\left|L_{i,j}\right|\le N_{max}$, where $M_L$ is the number of samples for each problem $q_i$;
• For each sample $L_{i,j}$, calculate the vector of ELA features $C_{i,j}$ and find the vector $P_{i,j}^{*}$ that is the best according to some metric $\varphi \left(L_{i,j},P\right)$ by solving the problem (3);
• Build a tuning model $\widehat{P}\left(C\right)$ using the set of known pairs $\left\{C_{i,j},P_{i,j}^{*}\right\}$, the total number of which is $M_Q\times M_L$.

Exploitation phase—at step 2a of the Bayesian optimization algorithm:

• Calculate the vector of features $C^r$ using the current training sample $L^r$;
• Predict the best hyperparameter values ${\widehat{P}}^r=\widehat{P}\left(C^r\right)$ using the tuning model built at step 4 of the exploration stage;
• Build the surrogate model ${\widehat{f}}^r\left(X\right)$ using the current training sample $L^r$ and the vector ${\widehat{P}}^r$ of the predicted hyperparameter values.

The efficiency of hyperparameter prediction is generally defined by the set $Q$, the feature vector $C$, and the metric $\varphi \left(L,P\right)$. The exploration phase involves most of the computational expenses of solving the problem (3), which makes it practical to apply even costly metrics, such as the metric ${\varphi }_{CV}\left(L,P\right)$. As illustrated in Figure 3, the test problems set $Q$ can be refined and extended permanently, even during the exploitation phase. New metrics are proposed in Section 3 to speed up the quality estimation of surrogate models for hyperparameter tuning and prediction. The efficiency of hyperparameter prediction with the metric ${\varphi }_{CV}\left(L,P\right)$ and with the proposed new metrics will be examined in Section 4. The scope of the article does not include formation of a representative set $Q$ of the test problems, or identification of the most suitable vector $C$ of ELA features.

3. Exploratory Landscape Validation

This section presents new metrics for estimating the efficiency of vector $P$ on a given sample $L$. The new metrics evaluate not only the approximation accuracy of the surrogate models, as the metric ${\varphi }_{CV}\left(L,P\right)$ does, but also certain properties of them that could be crucial for the optimization algorithm’s performance. Since the metrics are developed on the basis of an ELA algorithm, namely the variability map (VM) algorithm [12], we refer to them as exploratory landscape validation (ELV) metrics.

There are known metrics for estimating the quality of surrogate models that are not based on approximation accuracy and can, hence, be considered ELV metrics. For example, the ranking preservation (RP) metric [11], which we will refer to as ${\varphi }_{RP}\left(L,P\right)$, estimates the level of preservation of comparative relations between pairs of objective values approximated by the given surrogate model. According to the authors, the metric ${\varphi }_{RP}\left(L,P\right)$ is calculated using an independent test sample of points, that is not practical to collect in the context of a strict budget for objective function evaluations. At the same time, when using interpolating approximation algorithms, such as GP, the metric value calculated for the training sample will be close to the maximum possible value. To be able to use the metric ${\varphi }_{RP}\left(L,P\right)$ for hyperparameter tuning, an algorithm for generating an extended training sample $L_{ext}$ is introduced in this paper. The proposed metrics, which are based on extended samples $L_{ext}$, are referred to as ${\varphi }_{RP}\left(L_{ext},P\right)$ and ${\varphi }_{AD}\left(L_{ext},P\right)$. The metric ${\varphi }_{AD}\left(L_{ext},P\right)$ is unique in the sense that it estimates the quality of surrogate models by directly comparing the VMs of those models with the VM of the training sample.

Section 3.1 starts with an improved algorithm for building a VM, followed by a discussion of VM quality assessment. Next, in Section 3.2, an algorithm for constructing an extended variability map (EVM) from a VM is suggested to enhance the reliability of landscape validation by extending the training sample. Based on that, new EVM-based landscape validation metrics ${\varphi }_{RP}\left(L_{ext},P\right)$ and ${\varphi }_{AD}\left(L_{ext},P\right)$ are proposed and explained in Section 3.3. Thus, the proposed ELV approach consists of the process of building VMs, extending VMs, and calculating values of one of the suggested metrics based on an extended training sample when solving the problem (3).

3.1. Variability Map of an Objective Function

VMs were first proposed to estimate ELA features of the function $f\left(X\right)$ based on a given sample $L$ [12]. A VM is built by collecting a set of triples $T=\left\{t_j,j\in \left(1:\left|T\right|\right)\right\}$ from the sample $L$, where $\left|T\right|$ is the total number of triples. Each triple $t_j=\left(i_1^j,i_2^j,i_3^j\right)$ is composed of points $\left(X^{i_1},X^{i_2},X^{i_3}\right)$ that are neighboring in $X$ space. The triples for each point of sample $L$ are collected from all of its neighboring points. Two points are considered neighbors if the distance between them is less than the maximum distance between pairs of the closest points in sample $L$, taken with some correction factor. Using the corresponding objective values $\left(f^{i_1},f^{i_2},f^{i_3}\right)$ of the collected triples, the pair of increment values $\left({\delta }_1^j,{\delta }_2^j\right)$ is calculated. The values ${\delta }_1^j,{\delta }_2^j$ characterize increments of the objective function between pairs of points $\left(i_1,i_2\right)$ and $\left(i_2,i_3\right),$ respectively. The set of increment values $\left\{\left({\delta }_1^j,{\delta }_2^j\right),j\in \left(1:\left|T\right|\right)\right\}$ then forms the VM and can be visually represented as a cloud of points on the plane $0{\delta }_1{\delta }_2$, as shown in Figure 4.

In cases where sample $L$ is irregular, such as when point density varies a great deal, using a max–min distance estimate for triple collecting, as suggested in [12], may not be a sustainable strategy. We propose a new algorithm for VM building based on angular ranges, which works better with irregular samples as well.

The new algorithm consists of the steps listed below:

• Select a random point $X^{i_2}$ from sample $L$ and find the closest point $X^{i_3}$:

  $$i_3=\underset{i_3\in \left[1:\left|L\right|\right],i_3\ne i_2}{\arg \min }{d}_{i_2,i_3},$$

  (5)

  where $d_{i_2,i_3}=\left|\right|X^{i_3}-X^{i_2}\left|\right|$ is the Euclidean distance between the points;
• Find all the points $\left\{X^k,k\in \left[1:\left|L\right|\right],k\ne i_2, k\ne i_3\right\}$ that satisfy the conditions:

  $$\begin{gathered} d_{k,i_2}<d_{k,i_3}, \\ {d}_{k,i_2}<{\overline{d}}_{i_2}, \\ {\alpha }_{k,i_2,i_3}\ge \pi /2, \end{gathered}$$

  (6)

  where ${\overline{d}}_{i_2}$ is the mean distance from $i_2$-th point to all the other points, ${\alpha }_{k,i_2,i_3}$ is the angle formed by $X^k$, $X^{i_2}$, $X^{i_3}$ points in the $X$ space. The first condition is a quick check that reduces the number of points for which the angle needs to be calculated;
• For each angle range $\left[{\alpha }^{-},{\alpha }^{+}\right]$ from the set of ranges $\left\{\left[90,120\right],\left[120,150\right],\left[150,180\right]\right\}$, that is formed by splitting the range $\left(90,180\right)$ into three equal ranges, do the following:
  • find a point $X^{i_1}\in \left\{X^k\right\}$ that has the minimal distance $d_{i_1,i_2}$ in that angle range:

    $$\begin{gathered} X^{i_1}=\underset{X^{i_1}\in \left\{X^k\right\}}{\arg \min }{d}_{i_1,i_2}, \\ {\alpha }^{-}<{\alpha }_{i_1,i_2,i_3}\le {\alpha }^{+}; \end{gathered}$$

    (7)
  • extend the set of triples $T$ with a new one $t=\left(i_1,i_2,i_3\right)$;
  • increase the distances $d_{i_1,i_2}$ and $d_{i_2,i_3}$ by a factor of 2 so that the other points are considered if $i_2$ is randomly selected in the next iterations.
• If the maximum number of triples $\left|T\right|$ is not reached, move to step 1.

The presented algorithm has several tunable parameters. The set of angle ranges is formed by splitting the range of allowed angles $\left(90,180\right)$ into three equal ranges. It is recommended to use a lower bound of at least 90 degrees as triples with smaller angles may not be informative for landscape analysis purposes [12]. At the same time the angle formed by a triple of points in the $X$ space is limited to 180 degrees. The number of splits determines the maximum number of triples to be collected for each considered point. It is recommended to increase the number of splits for bigger dimensions $\left|X\right|$. Both the lower bound of the allowed angle range and the number of splits are tunable parameters of the algorithm. The maximum number of triples $\left|T\right|$ is another parameter of the algorithm that is usually set as the multiple of the total number of points $\left|L\right|$, e.g., by multiplying $\left|L\right|$ by the number of angle ranges. Increasing the distance between the points in step 3c helps to avoid selecting the same points for the next triples. It is recommended to multiply the distance between $X^{i_1}$ and $X^{i_2}$ at least by 1.5 and completely remove $X^{i_3}$ from the $X^{i_2}$ neighbors (e.g., by setting the corresponding distance value to infinity).

In Figure 5, the difference between the old and the new algorithm for collecting VM triples based on an irregular sample of points is illustrated. The training sample of size $\left|L\right|=30$ is represented by black dots in the $X$ space with dimension $\left|X\right|=2$. The points of the collected triples are connected by multicolor lines. It can be seen from the figure that the new algorithm provides better “coverage” of the $X$ space with the same number of triples. In Figure 5b, triples form connections between more distant points and fill the gaps in $X$ space caused by an irregular sample structure. The new algorithm generates triples with a lesser number of shared pairs of points, which is better for the landscape validation algorithm, further described below, and which is based on extended training samples.

The quality of a VM-building algorithm can be formally measured by the average distance between the points of the triples and average angles formed by those points. Formalizing the “coverage” quality and analyzing the connection between the quality of the VM and the efficiency of hyperparameter tuning or prediction is beyond the scope of this article.

3.2. Extended Variability Map of an Objective Function

A training sample-based landscape validation is not suitable for interpolation techniques like GP; at the same time, it is not practical to evaluate an expensive objective function for additional points in such cases. We propose the following method for extending the training sample $L$ using an extended version of a VM, built on that sample.

• Collect the set of triples $T=\left\{t_j,j\in \left(1:\left|T\right|\right)\right\}$ as it was described in Section 3.1;
• For each triple $t\in T$, where $t=\left(i_1,i_2,i_3\right)$, perform the following steps:
  • split the vector $X^{i_1}{X}^{i_2}$ into three parts by the points $X^{k_1}$ and $X^{k_2}$, so that the points $X^{i_1},X^{k_1}, X^{k_2},X^{i_2}$ are arranged in the given order on the same line in $X$ space, where $k_1,k_2>\left|L\right|$ for the new points $X^{k_1}, X^{k_2}$;
  • calculate the approximate objective values ${\tilde{f}}^{k_1}, {\tilde{f}}^{k_2}$ for the new points $X^{k_1}, X^{k_2}$, respectively, using a linear interpolation between the known points $\left\{\left(X^{i_1},f^{i_1}\right),\left(X^{i_2},f^{i_2}\right)\right\}$;
  • update the training sample $L_{ext}=L_{ext}{\displaystyle \cup }\left\{\left(X^{k_1},{\tilde{f}}^{k_1}\right),\left(X^{k_2},{\tilde{f}}^{k_2}\right)\right\}$;
  • update the extended set of triples $T_{ext}=T_{ext}{\displaystyle \cup }\left\{\left(i_1,k_1,k_2\right),\left(k_1,k_2,i_2\right)\right\}$;
  • repeat steps a-d for the vector $X^{i_2}{X}^{i_3}$.

The set of triples $T_{ext}=\left\{t_j,j\in \left(1:\left|T_{ext}\right|\right)\right\}$ and the corresponding increment values $\left\{\left({\delta }_1^j,{\delta }_2^j\right),j\in \left(1:\left|T_{ext}\right|\right)\right\}$ form the extended variability map (EVM). The EVM is based on an extended training sample $L_{ext}$, which includes the new points linearly interpolated between the points of triples $T$. The examples of the original training sample $L$ and the extended sample $L_{ext}$, built with the proposed algorithm, are shown in Figure 6. It is clear from Figure 6b that the new points are placed along the triples of the original sample presented in Figure 6a.

Note that the set of triples $T_{ext}$ does not include the original set $T$, while the sample $L_{ext}$ also includes the points from $L$. It is recommended to locate the new points $X^{k_1}$, $X^{k_2}$ closer to the points $X^{i_1},X^{i_2}$, e.g., by making logarithmic steps when splitting the vector $X^{i_1}{X}^{i_2}$, since it may positively affect the accuracy of landscape validation. The closer the new points are to the original ones, the stronger the requirements for the surrogate model to preserve comparative relations between pairs of those points.

In Figure 7, the VM of the original sample and the EMV of the extended sample can be seen. Since the new points are linearly interpolated, the additional points of the EVM (Figure 7b) are located on the diagonals ${\delta }_1={\delta }_2$ and ${\delta }_1=-{\delta }_2$.

3.3. Landscape Validation Metrics

To measure the landscape consistency between a surrogate model and the original sample $L^r$ using the extended training sample $L_{ext}^r$, the value of ${\varphi }_{RP}\left(L_{ext}^r,P\right)$ metric is calculated the following way:

• Build a surrogate model $\widehat{f}\left(X\right)$ using the training sample $L^r$ and the vector $P$ of hyperparameter values;
• Using the model $\widehat{f}\left(X\right)$, calculate ${\widehat{f}}^i$ values for all the points $X^i$of the extended sample $L_{ext}^r$, where $i\in \left(1:\left|L_{ext}\right|\right)$. The corresponding values form the sample ${\widehat{L}}_{ext}^r$;
• Given the samples $L_{ext}^r$ and ${\widehat{L}}_{ext}^r$, that have the common $X^i$ values, calculate the ranking preservation metric:

  $${\varphi }_{RP}\left(L_{ext}^r,P\right)=\frac{1}{\left(\begin{array}{c}\left|L_{ext}^r\right|\\ 2\end{array}\right)}{\displaystyle \sum }_{i=1}^{\left|L_{ext}^r\right|}{\displaystyle \sum }_{j=i+1}^{\left|L_{ext}^r\right|}\{\begin{array}{l}1, if comp\left(f^i,f^j\right)=comp\left({\widehat{f}}^i,{\widehat{f}}^j\right)\hfill \\ 0, otherwise\hfill \end{array},$$

  (8)

  where $comp\left(f^i,f^j\right)$ is the result of the comparison of $f^i$ and $f^j$ values with the possible outcomes: less, equal, more.

Note that $L_{ext}^r$ includes both the original values $f^i$ from L and the linearly interpolated values ${\tilde{f}}^i$. The value of the metric ${\varphi }_{RP}\left(L_{ext}^r,P\right)$ in the range $\left(0,1\right)$ measures the ratio of pairwise comparisons of objective values that are preserved by the model $\widehat{f}\left(X\right)$.

Figure 8 shows an example of models with different levels of ranking preservation. The training sample points are shown in black, while the predicted model values are shown in blue. The highest level of ranking preservation is achieved for Model 1, shown in Figure 8a, as for any pair of points $x^1$, $x^2$ in the given range the ranking preservation condition $comp\left(f^1,f^2\right)=comp\left({\widehat{f}}^1,{\widehat{f}}^2\right)$ is met. Model 2 in Figure 8b is violating the rankings near the edge points of the training sample. Although Model 1 preserves the ranking better, both Model 1 and Model 2 will have the highest value of the metric ${\varphi }_{RP}\left(L,P\right)$ when estimated based on the training sample $L$ only.

Figure 9 illustrates the idea of using the extended training sample to identify the level of ranking preservation. Although Model 2, shown in Figure 9b, has a better accuracy on the training sample, the violation of extra point ranking indicates that the model may not be suitable for optimization purposes. Although the metrics ${\varphi }_{RP}\left(L_{ext},P\right)$ and ${\varphi }_{CV}\left(L,P\right)$ are correlated, they are not identical, meaning that models with similar accuracy may preserve the landscape of the training sample differently.

We propose another landscape validation metric ${\varphi }_{AD}\left(L_{ext},P\right)$, called the angular divergence (AD) of the EVM. The metric is based on the extended sample $L_{ext}$, the extended set of triples $T_{ext}$ and the corresponding increment values $\left\{\left({\delta }_1^j,{\delta }_2^j\right),j\in \left(1:\left|T_{ext}\right|\right)\right\}$. Instead of directly measuring the consistency of the $\widehat{f}\left(X\right)$ model’s landscape using the original sample $L$, it measures the consistency of the corresponding EVMs in the following way:

• Build a surrogate model $\widehat{f}\left(X\right)$ using the training sample $L^r$ and the vector $P$ of hyperparameter values;
• Using the model $\widehat{f}\left(X\right)$, calculate ${\widehat{f}}^i$ values for all the points $X^i$of the extended sample $L_{ext}^r$, where $i\in \left(1:\left|L_{ext}^r\right|\right)$. The calculated values form the sample ${\widehat{L}}_{ext}^r$;
• For all the triples from $T_{ext}^r$, calculate the increment values $\left\{\left({\widehat{\delta }}_1^j,{\widehat{\delta }}_2^j\right),j\in \left(1:\left|T_{ext}^r\right|\right)\right\}$ using the sample ${\widehat{L}}_{ext}^r$;
• Assuming each pair of the increment values $\left({\delta }_1^j,{\delta }_2^j\right)$ correspond to the vector ${\mathsf{\Delta }}^j$ in the ${\delta }_1{\delta }_2$ space, and each pair of the values $\left({\widehat{\delta }}_1^j,{\widehat{\delta }}_2^j\right)$—to the vector ${\widehat{\mathsf{\Delta }}}^j$, calculate the angular divergence metric based on the cosine similarity of the vectors ${\mathsf{\Delta }}^j$ and ${\widehat{\mathsf{\Delta }}}^j$:

  $${\varphi }_{AD}\left(L_{ext}^r,P\right)=\frac{1}{\left|T_{ext}^r\right|}{\displaystyle \sum }_{j=1}^{T_{ext}^r}\frac{{\mathsf{\Delta }}^j·{\widehat{\mathsf{\Delta }}}^j}{||{\mathsf{\Delta }}^j||||{\widehat{\mathsf{\Delta }}}^j\left|\right|},$$

  (9)

  where ${\mathsf{\Delta }}^j·{\widehat{\mathsf{\Delta }}}^j$ is a dot product of the corresponding vectors.

The metric ${\varphi }_{AD}\left(L_{ext},P\right)$ is calculated as an average angle of rotation of EVM points around the central point of the $0{\delta }_1{\delta }_2$ plane. The rotation angle is determined by the model’s ability to preserve the ratio between the ${\delta }_1^j$and ${\delta }_2^j$ values, and not the absolute magnitude of those values. In Figure 10, the example of the AD calculation of a single triple for models with different levels of ranking preservation is shown. Model 1, presented in Figure 10a, preserves the signs of the increment values, so the angle between the EVM points that correspond to the sample and the model is relatively small. Model 2, shown in Figure 10b, alters the sign of the second increment value, hence the angular divergence for the corresponding EVM point is around 90 degrees. The signs of both increments are altered by Model 3 shown in Figure 10c, which results in an even larger angular divergence. Figure 11 illustrates the angular divergence for the whole EVMs built for GP models with different values of the Matern kernel parameter. The extended sample increments are shown in black points, while the increments calculated for the models are shown in blue. The value of the metric ${\varphi }_{AD}\left(L_{ext},P\right)$ of the first model shown in Figure 11a is lower than for the second model shown in Figure 11b (13 and 19 degrees on average, correspondingly).

4. Computational Experiment

In this section, the results of computational experiments performed with the proposed landscape validation metrics are presented. The baseline is set by estimating the Bayesian optimization efficiency with fixed hyperparameter values. In the first experiment, the efficiency of Bayesian optimization with hyperparameter tuning based on different landscape validation metrics is estimated. The second experiment evaluates the efficiency of the hyperparameter prediction approach based on ELA-ML approach with different landscape validation metrics.

4.1. General Setup

The Bayesian optimization efficiency is estimated on the set $Q=\left\{q_i,i\in \left(1:M_Q\right)\right\}$ of optimization problems that are based on the BBOB set of 24 test functions, that are widely used for optimization efficiency studies [20]. The BBOB set is accessed by using Python interface of the IOHexperimenter package [21]. Each function in the BBOB set is available for arbitrary dimensions $\left|X\right|$ and the numbers of so-called instances obtained by a random shift in $X$ space. The set $Q$ of test optimization problems is composed of BBOB test functions that have dimensions $\left|X\right|=2,4,8$ and the fixed instance number. Hence the total number of test problems is $M_Q=72$.

To solve test problems $Q$, an open Python implementation of the Bayesian optimization algorithm in the ByesOpt package is used [22]. The surrogate models are built using the Matern kernel, which has a tunable parameter $\nu$ with a recommended value of ${\nu }^{rec}=2.5$ specified in the package. The set of allowed $\nu$ values is fixed to $D_{\nu }=\left\{0.5, 1.5, 2.0, 2.5, 3.0,inf\right\}$ for hyperparameter tuning and prediction. By default, the package uses the LCB acquisition function with $\kappa =2.576$ to select the next point for objective evaluation [22].

The set of ELA features $C=\left(c_1,\dots c_{84}\right)$ is used to categorize the test problems for hyperparameter prediction. The vector $C$ includes 41 features evaluated by the pFlacco package [23] and 43 features based on VM [12]. Only ELA features that are based on a fixed sample of points are used since additional objective evaluations are not permissible with a fixed computational budget. The ELA features that require cell mapping of the search space are also excluded. Due to the exponential growth of the total number of cells with the dimension $\left|X\right|$, the sample sizes will not be sufficient for cell mapping feature calculation.

In the following sections, the method of solving the tuning problem (3) is referred to as a strategy of the Bayesian optimization algorithm. The efficiency of the following strategies is analyzed:

• Fixed vector $P$—no hyperparameter tuning;
• Hyperparameter tuning with the considered metrics:

  -

  metric ${\varphi }_{CV}\left(L,P\right)$ with 5 folds (see Equation (4) in Section 2.2);

  -

  metric ${\varphi }_{RP}\left(L_{ext},P\right)$ (see Equation (8) in Section 3.3);

  -

  metric ${\varphi }_{AD}\left(L_{ext},P\right)$ (see Equation (9) in Section 3.3).

• Predicted vector $\widehat{P}$ (see Section 2.3);
• The metrics ${\varphi }_{RP}\left(L_{ext},P\right)$ and ${\varphi }_{AD}\left(L_{ext},P\right)$ will be referred to as ${\varphi }_{RP}\left(L,P\right)$ and ${\varphi }_{AD}\left(L,P\right),$ respectively, to simplify the experiment description.

In the experiments, the efficiency of the Bayesian optimization algorithm with a given strategy is measured by the solution’s quality and by the computational cost of solving a test problem. The solution’s quality is determined by the best objective value found during the run. The computational cost is assessed by estimating the average time required to solve a test problem, which includes the time spent on building surrogate models and calculating values of the selected metric. Since test problems are used, the time spent on objective evaluation is negligible. The time required to build a tuning model is not taken into account when using the hyperparameter prediction strategy. The experiments were performed on a computer with an Intel E5-2643 processor.

4.2. Bayesian Optimization with Hyperparameter Tuning

To compare the efficiency of Bayesian optimization with fixed hyperparameter values and hyperparameter tuning, the experiment with the following steps is performed.

• For each problem $q_i\in Q$ with the objective function $f_i\left(X\right)$, generate the number of random initial samples $\left\{L_{i,j}^0, j\in \left(1:10\right)\right\}$, each sample of size $\left|L_{i,j}^0\right|=5\left|X\right|$;
• Using each initial sample $L_{i,j}^0$, perform iterations $r\in \left(1:10\left|X\right|\right)$ of the Bayesian optimization algorithm, as described in Section 2.1;
• Using each initial sample $L_{i,j}^0$, perform iterations $r\in \left(1:10\left|X\right|\right)$ of the Bayesian optimization algorithm with hyperparameter tuning, as described in Section 2.2, with each of the following metrics: ${\varphi }_{CV}\left(L_{i,j}^r,P\right), {\varphi }_{RP}\left(L_{i,j}^r,P\right)$ and ${\varphi }_{AD}\left(L_{i,j}^r,P\right)$. The best hyperparameter values are selected from the set $D_{\nu }$;
• Estimate the average of the best objective values ${\overline{f}}_i^{*}$ found in 10 random runs with fixed hyperparameter values in step 2 and with hyperparameter tuning based on the metrics ${\varphi }_{CV}\left(L_{i,j}^r,P\right), {\varphi }_{RP}\left(L_{i,j}^r,P\right)$ and ${\varphi }_{AD}\left(L_{i,j}^r,P\right)$ in step 3:

  $${\overline{f}}_i^{*}=\frac{1}{10}{\displaystyle \sum }_j{f}_{i,j}^{*},$$

  (10)

  where $f_{i,j}^{*}$ is the best objective value found for the problem $q_i$ in the $j$-th run;
• For each problem $q_i$ select the metric ${\varphi }_i\left(L,P\right)$ with which the best objective value was found on average.

In the described experiment, each of the 72 test optimization problems was solved using 10 random initial samples and four different strategies with the computational budget $15\left|X\right|$−2880 runs of the Bayesian optimization algorithm in total.

4.3. Bayesian Optimization with Hyperparameter Prediction

The experiment aims to measure the efficiency of using the hyperparameter prediction approach based on the ELA-ML framework by performing the following steps for each optimization problem $q_i$, $i\in \left(1:72\right)$. The experiment is split into exploration and exploitation phases as it is described in Section 2.3.

Exploration phase—evaluate ELA features, find the best hyperparameter values:

• Remove problems from the set $Q$ that are based on the same BBOB function as the current problem $q_i$, including those with different dimensions $\left|X\right|$, so that the remaining problems compose the set $Q_i=\left\{q_k,k\in \left(1:69\right)\right\}$;
• For each problem $q_k$ generate the random samples $\left\{L_{k,s},s\in \left(1:300\right)\right\}$ with different sizes $\left|L_{k,s}\right|\in \left(5\left|X\right|, 15\left|X\right|\right)$. Within the given range, 20 sample sizes are chosen and 15 random samples of each size are generated in $D_X$;
• For each sample $L_{k,s}$ calculate the vector of ELA features C_k,s, where $\left|C_{k,s}\right|=84$ and find the best hyperparameter values $P_{k,s}^{*}$ by solving the problem (3) with the set of allowed values $D_{\nu }$. As the hyperparameter efficiency metric use the metric ${\varphi }_k\left(L,P\right)$, which showed the best performance for the problem $q_k$ in Section 4.2;
• Use the set of pairs $\left\{C_{k,s},P_{k,s}^{*}\right\}$ as a training sample to build a tuning model ${\widehat{P}}_i\left(C\right)$ by using the random forest classifier implemented in the scikit-learn package [24].

Exploitation phase—use the tuning model ${\widehat{P}}_i$ for hyperparameter prediction:

• For the current problem $q_i$, generate the number of random initial samples $\left\{L_{i,j}^0, j\in \left(1:10\right)\right\}$, each sample of size $\left|L_{i,j}^0\right|=5\left|X\right|$;
• Using each initial sample $L_{i,j}^0$ perform iterations $r\in \left(1:10\left|X\right|\right)$ of the Bayesian optimization algorithm with hyperparameter prediction by using the tuning model ${\widehat{P}}_i$, as described for the exploitation phase in Section 2.3;
• Estimate the average of the best objective values ${\overline{f}}_i^{*}$ found in 10 random runs with hyperparameter prediction.

Each tuning model ${\widehat{P}}_i$ is built using 20,700 pairs of observations, all of which are collected for problems based on different BBOB functions. The experiment is structured in such a way that it is comparable to the cross-validation procedure. Each problem $q_i$ is excluded from the exploration phase to build a tuning model and is solved during the exploitation phase with the hyperparameter values predicted by that model. The hyperparameter values are predicted based on the ELA features similarity between the problem $q_i$ and the set of problems $\left\{q_k\right\}$ analyzed during the exploration phase. In our case, the accuracy of hyperparameter prediction for a particular problem $q_i$ depends, among other factors, on the presence of BBOB functions with a similar landscape in the set $\left\{q_k\right\}$.

4.4. Experimental Results

The results of the first experiment are summarized in

Table 1. The number of best-solved problems with the fixed hyperparameter value and the hyperparameter tuning approach using different metrics.

Bayesian Optimization Strategy	$\mathbf{Number}\mathbf{of}\mathbf{Problems}\mathbf{with}\mathbf{the}\mathbf{Best}{\overline{\mathit{f}}}_{\mathit{i}}^{*}$
$\mathrm{Fixed}\mathrm{vector}P$	12
Hyperparameter tuning with:
$\mathrm{metric}{\varphi }_{CV}$	24
$\mathrm{metric}{\varphi }_{RP}$	24
$\mathrm{metric}{\varphi }_{AD}$	21

. For each problem $q_i$, the best strategy is selected that provides the best value ${\overline{f}}_i^{*}$. For each considered strategy, the table shows the number of problems $q_i$ where the best result was found while using that strategy. It should be noted that for certain test problems, such as with the linear slope objective function, multiple strategies were able to find the optimal solution.

With the fixed hyperparameter values, the number of best-solved problems is the lowest, as expected. Using the proposed metrics for hyperparameter tuning leads to the best results on a wider range of problems. Based on the results, it can be assumed that different problems require different approaches to hyperparameter tuning, i.e., different metrics.

Table 2. The experimental results with the fixed hyperparameter value, the hyperparameter tuning approach using different metrics and the hyperparameter prediction approach.

Bayesian Optimization Strategy	$\mathbf{Number}\mathbf{of}\mathbf{Problems}\mathbf{with}\mathbf{the}\mathbf{Best}\mathbf{Value}{\overline{\mathit{f}}}_{\mathit{i}}^{*}$	$\mathbf{Average}\mathbf{of}\mathbf{Normed}\mathbf{Best}\mathbf{Values}{\overline{\mathit{f}}}_{\mathit{i}}^{*}$	Average Time for Solving a Problem, s
$\mathrm{Fixed}\mathrm{vector}P$	9	0.379	24
Hyperparameter tuning with:
$\mathrm{metric}{\varphi }_{CV}$	19	0.303	457
$\mathrm{metric}{\varphi }_{RP}$	13	0.305	211
$\mathrm{metric}{\varphi }_{AD}$	17	0.305	173
$\mathrm{Predicted}\mathrm{vector}\widehat{P}$	26	0.290	144

summarizes the results of the second experiment in the same manner. On top of that, the average of the best-found values ${\overline{f}}_i^{*}$ and the average time required to solve a test problem are provided. Since the scale of objective values of the test problems vary a great deal, the best-found values ${\overline{f}}_i^{*}$ were normed to the range $\left(0;1\right)$, so that 0 and 1 correspond to the best and worst values $f_{i,j}^{*}$ found for the problem $q_i$ in all the runs with different strategies.

It is evident that the proposed hyperparameter prediction approach provides the best results for a larger number of problems. The metric ${\varphi }_{CV}$ that requires cross-validation of the surrogate models is the most time-consuming but also the most accurate metric. However, with the proposed metrics ${\varphi }_{RP}$ and ${\varphi }_{AD}$ the results of comparable quality can be found with over 50% less effort. By using the hyperparameter prediction strategy, the problems can be solved with even less time (up to 70%), while the quality of the results is about 5% better than for the most accurate metric ${\varphi }_{CV}$.

$\mathrm{Fixed}$$\mathrm{metric}{\varphi }_{CV}$$\mathrm{metric}{\varphi }_{RP}$$\mathrm{metric}{\varphi }_{AD}$$\mathrm{Predicted}$When benchmarking optimization algorithms, average efficiency estimates are often not sufficient for making informed conclusions. Figure 12 presents the experimental results in the form of so-called performance profiles [25]. The performance profile of each strategy is constructed by calculating the number of test problems with a better value ${\overline{f}}_i^{*}\le {\overline{f}}^{\prime }$ for all possible values ${\overline{f}}^{\prime }\in \left(0,1\right)$. As a result, the performance profile plot shows the number of test problems as a function of the quality of the solution. For example, the two most effective strategies for finding near-best values ${\overline{f}}_i^{*}\le 0.1$ use the fixed vector $P$ and hyperparameter tuning with the metric ${\varphi }_{CV}$. Figure 13 presents the same performance profiles for all the algorithm runs with random initial samples (10 random runs for each of the 72 problems).

In Figure 14 the performance profile plots are shown for the number of test problems as a function of the time required to solve a problem. It can be clearly seen that the plot has «jumping» segments due to the time difference between solving the problems of different dimensions $\left|X\right|=2,4,8$. The computational complexity grows with the problem dimension slower for the metric ${\varphi }_{AD}$ than for the metric ${\varphi }_{RP}$. The hyperparameter prediction strategy is relatively expensive for smaller dimensions due to the computational costs involved in estimating the features vector $C$ and evaluating the tuning model $\widehat{P}\left(C\right)$. Using the fixed vector $P$ and selecting hyperparameter values with the metric ${\varphi }_{CV}$ are obviously the strategies with the best and the worst time efficiency, respectively. Figure 15 presents the same performance profiles for all the algorithm runs with random initial samples (10 random runs for each of the 72 problems).

5. Discussion

The presented exploratory landscape validation approach enables the finding pf effective hyperparameter values without relying on approximation accuracy estimations of surrogate models. In cases of limited computational budget and hence relatively small training samples, landscape validation metrics provide a faster way to estimate the efficiency of hyperparameter values. The presented metrics summarize the essential differences between the landscapes of surrogate models and training samples, that are generally characterized by ELA features. The study, however, has the following potential limitations. The experiments were performed with test optimization problems, but the applicability of the proposed metrics is mainly determined by the computational complexity of objective evaluation and hence the computational budget. In the event of a higher complexity of objective evaluation and relatively small budget, the computational impact even of the cross-validation metric may become negligible. On the other hand, if the budget is large enough, the increase in the optimization algorithm’s efficiency may not be sufficient to justify the computational cost of hyperparameter tuning with any of the considered metrics. Analyzing applicability conditions for different tuning strategies when solving practical optimization problems could be the focus of future research. It is also promising to explore other ways of estimating the level of landscape feature preservation by surrogate models based on the known ELA methods, for example, by direct comparison of ELA feature vectors. Landscape validation can be based on various ways of measuring the differences between the variability maps of samples and surrogate models.

It was shown that landscape validation metrics can be used to both find and predict the best hyperparameter values during Bayesian optimization. Another potential extension of this work is to analyze the efficiency of hyperparameter prediction with different ELA algorithms and approximation algorithms used for building a tuning model. The efficiency of hyperparameter prediction is also affected by the level of similarity between the problems being solved at exploration and exploitation phases, which is difficult to formalize. One possible approach would be to use the ELA features of the test problem set to define the region of allowed feature values for the exploitation-phase problems. If the new problem has unrelated feature values, then the efficiency of the predicted hyperparameter values cannot be guaranteed; therefore, the default hyperparameter values should be used for that problem. In such cases, the test problem set should be refined to keep it representative of the problems solved during the exploitation phase.

6. Conclusions

The article considers different approaches to improving the efficiency of Bayesian optimization algorithms by selecting the best hyperparameter values of the surrogate modeling algorithm. The optimal vector of hyperparameter values is found based on a hyperparameter efficiency metric, which defines the way of measuring the quality of a surrogate model built with different vectors. The hyperparameter tuning problem is being solved at each iteration of Bayesian optimization, so using computationally demanding metrics may lead to a significant increase in the time spent solving the problem.

When solving computationally expensive optimization problems, the number of objective evaluations allowed is relatively small, as well as the size of the training sample for building a surrogate model. The commonly used efficiency metric for such cases is the cross-validation score, which requires building multiple surrogate models on different subsets of the training sample. In this article, a new approach is introduced called exploratory landscape validation (ELV), which includes the proposed hyperparameter efficiency metrics to assess the quality of surrogate models without considering the approximation error estimates. The experiments showed that hyperparameter tuning with the new metrics can provide solutions of comparable quality with less than half the time required when using the cross-validation metric. The experimental results also indicate that different metrics provide the best solutions for different optimization problems.

Another way of reducing the costs of hyperparameter tuning is to build a model that predicts the best hyperparameter values during Bayesian optimization based on ELA features estimated from the training sample. The hyperparameter prediction approach is based on collecting a test set of problems, estimating the ELA features of those problems, finding the best hyperparameter values according to an efficiency metric, and building a tuning model. In general, each optimization problem has its own best-suited efficiency metric for hyperparameter tuning. In the computational experiment, the tuning models are built to predict the hyperparameter values that are most effective according to the metrics chosen individually for each test optimization problem. With the suggested hyperparameter prediction approach and individual efficiency metrics, better-quality solutions were found in less than 70% of the time needed by a hyperparameter tuning approach based on cross-validation score. Even though additional computational expenses are required to create a tuning model, they are insignificant when compared to potential permanent improvement in the optimization algorithm’s efficiency.