A black-box optimization method with polynomial-based kernels and quadratic-optimization annealing

Here, we describe the details of the surrogate model construction for the black-box function (exploitation in a Bayesian optimization context). As noted above, the surrogate model constructed here must be a second-order polynomial to utilize QUBO solvers. Using a relatively low-order polynomial also contributes to circumventing the over-fitting of the constructed model parameters, as the training data is insufficient, which is typical for most BBO problems.

Now, let our first kernel function be:

Note that the inverse operation involved in Eq. (7) seems expensive when the number of samples $n$ becomes larger. In the present study, $(\mathbf{K}^{*}_{\mu}+\lambda\mathbf{I})^{-1}$ is sequentially updated by using the Woodbury formula [17] each time a new input-output pair is added at the end of optimization cycle (4) in Fig. 1(a). This suppresses the computational cost by a factor of $\mathcal{O}(n^{2})$, and the overall cost of computing Eq. (7) is $\mathcal{O}(n)$.

Finally, the coefficient matrix can be obtained as $\hat{\mathbf{Q}}_{\mu}=\sum_{i=1}^{n}\hat{c}_{i}\mathbf{x}_{i}^{*}\mathbf{x}_{i}^{*T}$, $\hat{\mathbf{q}}_{\mu}=\sum_{i=1}^{n}\hat{c}_{i}\mathbf{x}_{i}^{*}$, and the final surrogate model $\hat{f}_{\mu}$ for the black-box function $f$ based on the given dataset is:

We can remove the constant $c_{\mu}$ as this does not generally affect optimization.

Both $\hat{f}_{\mu}(\mathbf{x})$ in Eq. (8) and $\hat{f}_{\sigma}(\mathbf{x})$ in Eq. (12) are second-order polynomials, thereby the acquisition function in Eq. (9) with non-zero $\beta$ being QUBO-compatible. Also, all the coefficients in the polynomial can be analytically obtained rather than iteratively fitted, and this point is advantageous when constructing a surrogate model with a relatively large input dimension.

Kernel-QA inherently requires the variables to be binary to utilise a fast QUBO solver for the model function optimization. However, with appropriate variable encoding, real and integer variables as well as their mix can be considered. Conversion from non-binary to binary variables can be performed either before or after the construction of the model function $\hat{f}(\mathbf{x})$, which are named respectively a priori and a posteriori conversions here.

In the present study, a priori conversion is chosen due to its straightforwardness, with a domain-wall encoding method. A real variable $x$ with lower and upper bounds $(x_{lower},x_{upper})$ is discretized onto a vector $\mathbf{x}_{disc}$ of a size $n_{bins}$. The spacing of bins $\Delta$ can be either uniform or non-uniform. For example, a variable $x$ with $(x_{lower},x_{upper})=(0.0,1.0)$ is discretized into $n_{bins}=5$ bins, and $\mathbf{x}_{disc}=[x_{disc,1},\cdots,x_{disc,n_{bins}}]=[0.0,0.25,0.50,0.75,1.0]$.

In the a priori variable conversion method, variable conversion (encoding) from real/integer values to binary vectors is performed before the construction of the surrogate model at (1) in Fig. 1(a). Thus, any inputs to the constructed surrogate model in Eq. (9) are binary values (each non-binary value is converted to a binary vector with the size $n_{bins}-1$). Also, variable conversion (decoding) from binary vectors to real/integer values is required before evaluating the black-box function at (3) in Fig. 1(a). Evaluation of other encoding methods, such as one-hot encoding and binary encoding with/without a posteriori conversion, is an interesting topic, which is out of scope in this study.

In kernel-QA with $\beta=0$, the construction of an appropriate surrogate model (e.g., Eqs. (8)) is the key. Here, the model does not need to be highly accurate: With the complex black-box function and relatively small dataset, such a highly accurate model cannot be obtained, especially for large-dimension problems. In the BBO context, an appropriate model yields a positive correlation with the black-box function. Such correlation characteristics of the model may be assessed by the cross-correlation between the true output values of the black-box function in the training dataset and the values predicted by the trained surrogate model. Typically, such a correlation coefficient should be positive, preferably greater than 0.5 throughout the cycles.

Sometimes, the output values of the black-box function $f(\mathbf{x})$ may have an extensive dynamic range, and the output is sensitive to the input values. Such situations can be identified by looking at the initial data set $(\mathbf{x}^{*},\mathbf{y}^{*})$. Constructing an “appropriate” model for such $f(\mathbf{x})$ can be difficult and often requires some feature scaling on a par with the usual machine learning, for example.

Additional care must be taken for BBO purposes. Minimization problems focus on the model’s behavior at relatively small output values. Thus, the model must positively correlate with the black-box function among samples with relatively small output values even though other samples hold large output values. To construct an effective surrogate model for the above-mentioned situation, we consider a transformation that diminishes returns for large values of $f(\mathbf{x})$ and encourage the optimizer to explore regions of the input space where $f(\mathbf{x})$ is smaller. Such transformation is beneficial if there are local minima in the original $f(\mathbf{x})$ that you would like to avoid or de-emphasize in favor of a deeper global search.

The exponential transformation in Eq. (16) magnifies small output values of $f(\mathbf{x})$ while compressing large ones. Such transformation would be helpful in problems where the objective function has a long tail, or the optimizer needs to focus on further reducing already small values of $f(\mathbf{x})$.

In the present study, we consider the following artificial landscapes to assess the proposed method, kernel-QA. These functions are for the input of arbitrary dimensions $d$ and are formulated as follows.

• •

  The Rosenbrock function:

  $$f(\mathbf{x})=\sum_{i=1}^{n-1}\left[(1-x_{i})^{2}+100(x_{i+1}-x_{i}^{2})^{2}\right],$$

  (18)

  Global minimum: $f(1,1,\cdots,1)=0$.

• •

  The Rastrigin functions:

  $$f(\mathbf{x})=10n+\sum_{i=1}^{n}\left[x_{i}^{2}-10\cos(2\pi x_{i})\right],$$

  (19)

  Global minimum: $f(0,0,\cdots,0)=0$.

Depending on the purpose of the assessment, the input dimensions of $d=5$–$80$ are considered for real-variable cases, and $d=40$–$640$ are considered for binary-variable cases. Note that for the kernel-QA, due to the conversion from $d$ real variables (see Sec. 2.3), the variables $\mathbf{x}$ considered in the surrogate model (Eq. (9)) are binary. The number of such binary variables that the surrogate model and its optimization need to deal with is $d_{B}=d(n_{bins}-1)$ in the present study, where all the real variables are encoded using the $n_{bins}$ uniform bins. For the binary-variable cases, such variable conversion is unnecessary.

The kernel-QA is parametrically assessed in the present study for the artificial landscapes described in Sec. 3.1. For reference, results from typical Bayesian optimization are also compared. Here, optimization conditions considered in the assessments are described.

A summary of optimization conditions considered in the present assessments is shown in Table 1. For both kernel-QA and Bayesian optimization, each artificial landscape in Fig. 2 (with the shown input dimensions, i.e., the number of real/binary variables, $d$) is regarded as black-box and optimized by considering each condition shown in Table 1. Also, identical initial exploration data are used for kernel-QA and Bayesian optimization cases, whose length (the number of input-output pairs to/from the black-box function) is denoted by $n_{init}$. The initial exploration data are constructed based on randomly chosen input vectors. For each BBO method, optimization runs are repeated ten times independently with a different initial seed for the initial exploration data to collect some statistics.

For kernel-QA, the same annealing timeout of 5 seconds is considered for all conditions, and Fixstars Amplify Annealing Engine (Amplify AE) [12] is used for the optimization of the acquisition function. Here, Amplify AE is a GPU-based Ising machine that is available to the public. It can accept QUBO problems with up to 256,000 bits. Thus, many typical BBO problems would be solvable regarding decision variable dimensions. Also, for the variable conversion required for the real decision variable cases, the a priori conversion (Sec. 2.3) is applied. These conditions are also applied to the present FMQA runs considered for some assessments. As for the model function for kernel-QA, Eq. (9) is considered with $\beta=0$ unless otherwise noted, without including the linear terms ($\gamma=0$ in Eqs. (3) and (11)) for simplicity. For all the cases, the regularization parameter $\lambda=1$.

As for Bayesian optimization, the expected improvement (EI) is used for the acquisition function whose coefficients are determined analytically [19]. Also, the minimization of the acquisition function is performed by using the limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) algorithm [20], which would cost $\mathcal{O}(dn^{2}T)$, where $T$ is the number of iterations required. These settings are typical in Bayesian optimization and implemented using the GPyOpt library [21] in the present study. Note that of all the parameters shown in Table 1, only $d$, $(x_{low},x_{up})$ and $n_{init}$ are directly relevant to Bayesian optimization.

Ten independent runs (different initial data samples) are performed for each optimization case, and their average and standard deviation are discussed in the present assessment. Note that identical initial training data samples are used between different optimization methods.

Kernel-QA and Bayesian optimization are performed to find an optimal input for “black-box” functions, the Rosenbrock and Rastrigin functions, with the real-number input dimension $d=5$ under the condition labeled as r5h in Table 1. The evolution of the objective function values with the found best solution at $i$-th cycle, $f(\mathbf{x}_{best,i})$, averaged over ten independent optimization runs, is shown in Fig. 3. Generally, all optimization methods show reasonable trends. That said, the results clarify the different characteristics of different optimization methods.

For the Rosenbrock function (no local minima), after the 1000-th cycle, the values of $f(\mathbf{x}_{best,1000})$ averaged over the ten runs are 0.7 for kernel-QA, 4.9 for Bayesian optimization, and 13.4 for FMQA. Clearly, FMQA shows a digit worse result, and this is potentially because the stochastic gradient descent used in FM with the relatively small dataset might have resulted in local optimal model coefficients as discussed in Eq. (7), and such behavior does not seem beneficial for the Rosenbrock’s landscape, where the global minimum lies in a very narrow valley.

As for the Rastrigin function in Fig. 3(b), Bayesian optimization yields a rapid decrease of $f(\mathbf{x}_{best,i})$ at the early phase of optimization cycles (say first 100 cycles), but the slope of evolution becomes less than the other methods. Kernel-QA and FMQA catch up with Bayesian optimization after approximately 100 cycles and continue to find better solutions that yield smaller $f(\mathbf{x})$. The values of $f(\mathbf{x}_{best,1000})$ averaged over the ten runs are 1.4 for kernel-QA, 3.1 for Bayesian optimization, and 1.3 for FMQA.

An interesting observation is shown in Fig. 4 which plots the distance between the found best solution at $i$-th cycle $\mathbf{x}_{i,best}$ and the true solution $\mathbf{x}_{true}$, computed as $\|\mathbf{x}_{i,best}-\mathbf{x}_{true}\|$ for the Rastrigin case. The evolution shows that the solutions found by kernel-QA are actually closer to the true solution than the ones found by Bayesian optimization as soon as the optimization cycles are started. On the other hand, in terms of the $f(\mathbf{x}_{best,i})$ shown in Fig. 3(b), Bayesian optimization shows a rapid decrease and seems outperforms kernel-QA for the first 100 cycles. This conflicting observation reveals kernel-QA and FMQA successfully avoid being trapped in local minima for the Rastrigin function. Such a tendency to avoid local minima in kernel-QA is considered due to the low-order polynomial in its surrogate model, which would circumvent the overfitting of the model parameters with a relatively small dataset typical in BBO problems. A similar trend is shown in Figs. 3(b) and 4 for FMQA, which also utilizes a second-order polynomial in the surrogate model, ensures this insight.

Sec. 2.4 described the exponential transformation of the training data ($f(\mathbf{x})$ output values) to help the surrogate model be constructed appropriately for optimization even when the black-box function yields an extensive dynamic range. The exponential transformation requires a transformation parameter $c_{m}$, which can be determined based on the initial training data, typically, the averaged output value $\alpha_{exp}\langle\mathbf{y}_{init}\rangle$ of the samples in the initial training dataset with $\alpha_{exp}$ being unity by default. This section discusses the effect of this transformation and the sensitivity of $c_{m}$ (or $\alpha_{exp}$) on the overall optimization performance.

Figure 7 shows the optimization history with kernel-QA under the conditions r80 and b640 for the Rosenbrock and Rastrigin functions. Four $\alpha_{exp}=0.1,0.5,1.0,2.0$ and one condition without the exponential transformation are considered for each case. The results with $\alpha_{exp}=1.0$ are identical to the ones for r80 and b640 in Figs. 5 and 6. There are substantial improvements for all the cases using this transformation method with basically any $\alpha_{exp}$. The unity $\alpha_{exp}$ (default in the present study) seems quite a decent choice. However, the results show that twice or half of the value also yields reasonable (or better) optimization performance.

Figure 1(a) shows that serial optimization methods require initial training data to construct the first model function. A relatively sizeable initial data size may result in a “better” model function for the first few optimization cycles. In contrast, a smaller initial data size reduces the cost of evaluating black-box functions during the initial data construction. Also, with a smaller initial data size, the ratio of newly added data samples during the optimization cycles to the total data becomes more prominent. This may result in faster convergence as, at the same cycles, the surrogate model fits better to the samples obtained during the optimization cycles rather than the initial data, for example, obtained randomly. The number of samples in the initial training data in the present study is set to be $b_{init}=10$ as summarized in Table 1. Here, the effect of initial data size on the optimization progress is discussed.

Figure 8 shows optimization history with different initial data sizes $n_{init}=2$, $10$ and $100$. While Fig. 8 shows the optimization results for the Rastrigin function, the trend for the corresponding Rosenbrock function is very similar. As clearly shown, the choice of $n_{init}$ does not substantially affect the overall optimization performance. That said, the result with $n_{init}=100$ shows lagged evolution of $f(\mathbf{x}_{best,i})$ reduction for the binary-variable problem in Fig. 8(b). Also, such $n_{init}$ requires more black-box evaluations, which is impractical. While $n_{init}=10$ is considered the default in the present experiments, the optimization cycles in kernel-QA can start with as small as $n_{init=2}$ samples without affecting optimization performance.

The assessments for kernel-QA hereinbefore utilized a surrogate model (Eq. (8)), rather than an acquisition function (Eq. (9)). However, for some BBO problems, it may be beneficial to consider both Eqs. (8) and (12) in the context of the acquisition function. The balance of the expected value (exploitation) and uncertainty (exploration) is controlled by $\beta$ in Eq. (9). when $\beta=0$, Eq. (9) is the surrogate model.

Figure 9 compares evolution of $f(\mathbf{x}_{best,i})$ for the Rosenbrock and Rastrigin functions with zero and non-zero $\beta$ for the conditions r10, r80, b80 and b640 summarized in Table 1. The considered non-zero $\beta$ values are 0.0001, 0.001, and 0.01, where $\beta=0.01$ means more emphasis on exploration than $\beta=0$ (default value). The effect of $\beta$ seems more prominent, in both positive and negative manner, for the larger dimension, r80 than r10 for the real-variable case, and b640 than b80 for the binary-variable case.

For the more significant dimension cases (r80 and b640), the optimization performance with $\beta=0.01$ (and larger $\beta$, expectedly) seems poorer regardless of functions or variable types. This result suggests the adverse effect of (excess) exploration for larger dimension problems, where the number of samples needed to sufficiently “cover” the search space grows rapidly with dimensionality. As for smaller non-zero $\beta$, its influence on overall optimization performance differs depending on the functions and variable types, and the trend does not seem consistent entirely. Such inconsistency is due to approximating the standard deviation in Eq. (12). However, this inconsistency reduces towards the end of optimization cycles, and at the 1000-th cycle, the values of $f(\mathbf{x})$ converge to similar values between $\beta=0,0.0001$, and $0.001$ for each case.

Although the present experiment does not demonstrate the consistent effectiveness of $\beta$, the result shows that the optimization performance could be further improved at optimization cycles $i\ll 1000$ for relatively large dimension problems if an appropriate $\beta$ value can be used. Alternatively, the QUBO-compatible formulation of $\sigma$ (simplification from Eq. (12) to Eq. (10)) could be explored to achieve consistent performance in future studies.