StepDIRECT - A Derivative-Free Optimization Method for Stepwise Functions

A common approach for the decision-making problem based on data-driven tools is to build a pipeline from historical data, to predictive model, to decisions. A two-stage solution is often used, where the prediction and the optimization are carried out in a separate manner [5]. Firstly, a machine learning model is trained to learn the underlying relationship between the controllable variables and the output. Secondly, the trained model is embedded in the downstream optimization to produce a decision. We assume that the regression model estimated from the data is a good representation of the complex relationship. In this paper, we consider the problem of optimizing a tree ensemble model such as random forests and boosting trees, where the predictive model has been learned from historical data.

The tree ensemble model combines predictions from multiple decision trees. A decision tree uses a tree-like structure to predict the outcome for an input feature vector ${\bf{x}}\in\mathbb{R}^{p}$. The $t$-th regression tree in the ensemble model has the following form

where $\Omega_{t,1},\ldots,\Omega_{t,M_{t}}$ represent a partition of feature space. We can see that $f_{t}({\bf{x}})$ is a stepwise function. The tree ensemble regression function outputs predictions by taking the weighted sum of multiple decision trees as

where $w_{t}$ is the weight for the decision tree $f_{t}({\bf{x}})$. We assume that parameters $c_{t,i},\Omega_{t,i}$ and $w_{t}$ have been learned from data, and we use StepDIRECT to optimize $f({\bf{x}})$. As we can see that $f_{t}({\bf{x}})$ is constant over each sub-region $R_{t,i}$; hence $f({\bf{x}})$ is a stepwise function. Formally, we have the following theorem.

We note that for this type of regression functions, we might get additional information about the objective function such as variable importance for input features of random forests [1].

In machine learning models, we often need to choose a number of hyper-parameters to get a high prediction accuracy for unseen samples [13]. For example, in training an $\ell_{1}$-regularized logistic regression classifier, we tune the sparsity parameter [19]. A common practice is to split the entire dataset into three subsets: training data, validation data, and test data [11]. For a given set of hyper-parameters, we train the model on the training data, then evaluate the model performance on the validation data. Some widely-used performance metrics include accuracy, precision, recall, $F_{1}$-score, and AUC (Area Under the Curve). The goal is to tune hyper-parameters for the model to maximize the performance on the validation data. The task can be formulated as a black-box optimization problem.

First, we start with a binary classifier. Suppose that we are given the $M$ training samples $\{({\bf{u}}_{1},{\bf{v}}_{1}),\ldots,({\bf{u}}_{M},{\bf{v}}_{M})\}$ and $N$ validation samples $\{({\bf{u}}_{1},{\bf{v}}_{1}),$ $\ldots,({\bf{u}}_{N},{\bf{v}}_{N})\}$ where ${\bf{u}}_{i}\in\mathbb{R}^{d}$ and ${\bf{v}}_{i}\in\{\pm 1\}$. For a fixed set of model parameters $\bm{\lambda}\in\mathbb{R}^{p}$, a classification model $h(\cdot;\bm{\lambda}):\mathbb{R}^{d}\to\{\pm 1\}$ has been learned based on the training data. When the performance measure is accuracy, the HPT problem is to determine $\bm{\lambda}$ that maximizes the test accuracy

where $({\bf{u}}_{i},{\bf{v}}_{i})$ is from the validation data. We can see that $F_{acc}(\bm{\lambda})\in\{0,\frac{1}{N},\ldots,\frac{N-1}{N},1\}$ for any $\bm{\lambda}$; hence we expect that the landscape of this function should have stepwise behavior. As an example, we plot in Figure 1 the landscape for the HPT logistic regression by tuning the $\ell_{1}$ parameter (denoted by $C$ in the $x$-axis) for training the a1a dataset [2].

The target function for other metrics can also take only a finite number of values. The observation still holds true for a multi-label classifier.

We follow the first step of the original DIRECT for initialization [8]. In this step, Algorithm 1 starts with finding $f_{1}=f({\bf{c}}_{1})$ and divides the hyper-cube ${\Omega}$ by evaluating the function values at $2p$ points ${\bf{c}}_{1}\pm\delta{\bf{e}}_{i},i\in\{1,\ldots,p\}$, where $\delta=\frac{1}{3}$ and ${\bf{e}}_{i}$ is the $i$-th unit vector. The idea of DIRECT is to select a hyper-rectangle with a small function value in a large search space; hence let us define

and the dimension with the smallest $s_{i}$ is partitioned into thirds. By doing so, ${\bf{c}}_{1}\pm\delta{\bf{e}}_{i}$ are the center of the newly generated hyper-rectangles. We initialize the sets ${\cal{P}}_{1},{\cal{I}}_{1},{\cal{C}}_{1}$, values $f_{{\cal{H}}_{i}}$ for every ${\cal{H}}_{i}\in{\cal{P}}_{1}$, and update $f_{min}=\min\{f_{{\cal{H}}_{i}}:{\cal{H}}_{i}\in{\cal{P}}_{1}\}$ and corresponding ${\bf{x}}_{min}$.

In this subsection, we propose a new criterion for StepDIRECT to select the next potentially optimal hyper-rectangles, which should be divided in this iteration. In the original DIRECT algorithm, every hyper-rectangle ${\cal{H}}_{i}$ is represented by a pair $(f_{{\cal{H}}_{i}},d_{i})$, where $f_{{\cal{H}}_{i}}$ is the function value estimated at the centre ${\bf{c}}_{i}$ of ${\cal{H}}_{i}$ and $d_{i}$ is the distance from the center of hyper-rectangle ${\cal{H}}_{i}$ to its vertices. The criterion to select hyper-rectangles, i.e., potentially optimal hyper-rectangles, for further divided is based on a score computed from $(f_{{\cal{H}}_{i}},d_{i})$. A pure local strategy would select the hyper-rectangle with the smallest value for $f_{{\cal{H}}_{i}}$, while a pure global search strategy would choose one of the hyper-rectangles with the biggest value $d_{i}$. The main idea of the DIRECT algorithm is to balance between the local and global search, which can be achieved by using a score weighting the two search strategies: $f_{{\cal{H}}_{i}}-Kd_{i}$ for some $K>0$. The original definition for potentially optimal hyper-rectangles ${\cal{H}}_{j}$ for DIRECT [16] is based on two conditions:

for some $\epsilon>0$.

An example for identifying these potentially optimal hyper-rectangles by the original DIRECT is given in Figure 2. Each dot in a two dimensional space represents a hyper-rectangle, three red dots with the smallest value for $f({\bf{c}}_{j})-Kd_{j}$ for each $K$ and a significant improvement (i.e., $f({\bf{c}}_{j})-Kd_{j}>f_{min}-\epsilon$) are considered as potentially optimal.

The proposed StepDIRECT searches locally and globally by dividing all hyper-rectangles that meet the criteria in Definition 2. A hyper-rectangle ${\cal{H}}_{i}$ is now represented by a triple $(f_{{\cal{H}}_{i}},d_{i},\sigma_{i})$, in which we introduce new notations for $f_{{\cal{H}}_{i}}$ and $\sigma_{i}$. We use a higher quality value $f_{{\cal{H}}_{i}}$ returned by a local search as a representative for each hyper-rectangle, and a flatness measurement $\sigma_{i}$ in the proposed criterion.

In Eq. (3.7), we introduce a new notation $\sigma_{j}$ to measure the local landscape of the stepwise function. Furthermore, we replace $f({\bf{c}}_{j})$ as in the original DIRECT by $f_{{\cal{H}}_{j}}$, which is computed with the help of local search. The value $f_{{\cal{H}}_{j}}$ for ${\cal{H}}_{j}$ generated during Steps a and c is better estimating the global solution of $f$ over ${\cal{H}}_{j}$ because $f_{{\cal{H}}_{j}}\leq f({\bf{c}}_{j})$. In Eq. (3.8), in order to eliminate the sensitivity to scaling issue, $|f_{min}|$ is replaced by the difference $|f_{min}-f_{median}|$ as suggested in [7].

Now we discuss how to compute $\sigma_{j}$. At the $k$-th iteration of the StepDIRECT algorithm, for hyper-rectangle ${\cal{H}}_{j}\in{\cal{P}}_{k}$ with center ${\bf{c}}_{j}$ and diameter $2d_{j}$, we define its set of neighborhood hyper-rectangles as

for some $\lambda>0$. Then, ${\cal{N}}_{j}$ can be further divided into two disjoint subsets ${\cal{N}}_{j}^{D}$ and ${\cal{N}}_{j}^{E}$ such that

The local variability estimator for hyper-rectangle ${\cal{H}}_{j}$ is defined as

where $0<\epsilon_{\sigma}<1$ is a small positive number to prevent $\sigma_{j}=0$ when ${\cal{N}}_{j}^{D}=\emptyset$. In our experiments, we set $\lambda=2$ and $\epsilon_{\sigma}=10^{-8}$ as default. The meaning of $\sigma_{j}$ can be interpreted as follows. For each hyper-rectangle ${\cal{H}}_{j}$, a large value for $\sigma_{j}$ indicates that it is likely to have more different function values in the neighbourhood of the center ${\bf{c}}_{j}$, which requires more local exploration. We propose to include the hyper-rectangle to the set of potentially optimal hyper-rectangles.

By Definition 2, we can efficiently find potentially optimal hyper-rectangles based on the following lemma.

We note that all proofs are delegated to the supplementary material. Definition 2 differs from the potentially optimal hyper-rectangle definition proposed in [16] (i.e., (3.5),(3.6)) in the following three aspects:

• •

  (3.7) and (3.8) include the local variability of $f$ on the hyper-rectangle. For the stepwise function, this quantity helps balance the local exploration and global search.

• •

  (3.8) uses $|f_{min}-f_{median}|$ to remove sensitivity to the linear and additive scaling and improve clustering of sample points near optimal solutions.

• •

  $f_{{\cal{H}}_{i}}$ can be different from $f({\bf{c}}_{i})$, i.e., $f_{{\cal{H}}_{i}}\leq f({\bf{c}}_{i})$.

Once a hyper-rectangle has been identified as potentially optimal, StepDIRECT divides this hyper-rectangle into smaller hyper-rectangles. We will take into account the side length and the variable importance measure of each dimensions.

In tree ensemble regression optimization, we can obtain the variable importance which indicate the relative importance of different features [1]. In general, the variable importance relates to a significant function value change if we move along this direction, and also the number of discontinuous points along the coordinate. As a result, we tend to make more splits along the direction with higher variable importance. Formally, let ${\bf{w}}=(w_{1},\ldots,w_{p})\in\mathbb{R}_{+}^{p}$ be the normalized variable importance with $||{\bf{w}}||_{1}=1$ and the length of the hyper-rectangle be ${\bf{l}}=(l_{1},\ldots,l_{p})\in\mathbb{R}_{+}^{p}$, then we define

as the relative importance of all coordinates for the selected hyper-rectangle. We choose the coordinate with the highest value $v_{i}$ to divide the hyper-rectangle into thirds.

If no variable importance measure is provided, then we take ${\bf{w}}=(1/p,\ldots,$ $1/p)\in\mathbb{R}_{+}^{p}$. The dividing procedure is the same as the original DIRECT by splitting the dimensions with the largest side length.

In this subsection, we introduce a randomized local search algorithm designed for minimizing stepwise function $f({\bf{x}})$ over a bounded box

It has been known that a DIRECT-type algorithm has a good ability to locate promising regions of the feasible space, and a good local search procedure can help to quickly converge to the global optimum [21, 20]. Since the function values over the neighborhood of a point are almost surely constant for a stepwise function, existing local searches based on a local approximation of the objective function might fail to move to a better feasible solution. Hence, we propose a novel local search method, shown in Algorithm 2. Different from the classical trust-region methods, Algorithm 2 will increase the stepsize when no better solution is discovered in the current iteration, i.e., $\tau>1$. This change is motivated by the stepwise landscape of the objective function.

We provide two different options to generate the search directions $\mathcal{D}$:

• (a)

  By coordinate strategy, for each axis $i$, we take two directions ${\bf{e}}_{i}$ and $-{\bf{e}}_{i}$ with a probability $w_{i}$ calculated from variable importance.

• (b)

  By random sampling from the unit sphere in $\mathbb{R}^{p}$.

The first option is more suitable when the discontinuous boundaries are in parallel to the axes, for example, ensemble tree regression functions. The second strategy works for general stepwise functions when the boundaries for each level set is not axis-parallel. We store the feasible sampled points ${\bf{x}}+\delta{\bf{d}}$ with their function values generated during local search in order to update the best function value $f_{{\cal{H}}_{i}}$ for each hyper-rectangle ${\cal{H}}_{i}$ for a new partition.

Now, we are able to provide the global convergence result for StepDIRECT. We assume that ${\bf{w}}=(1/p,\ldots,1/p)\in\mathbb{R}^{p}$ for simplicity and the following theorem is still valid as long as $w_{j}>0$ for every $j\in\{1,\ldots,p\}$.

We consider the minimization for Random Forest regression function over a bounded box constraint. We used the boston, diabetes, mpg and bodyfat data sets from the UCI Machine Learning Repository [6]. Table 1 provides the details of these four data sets.

We train the random forest regression function on these data sets with 100 trees and use default settings for other parameters in scikit-learn package [25]. For comparison, we run the following optimization algorithms: DIRECT [16]¹¹1https://scipydirect.readthedocs.io/en/latest/, DE [28]²²2https://docs.scipy.org/doc/scipy/reference/index.html, PSO [17]³³3https://pypi.org/project/pyswarms/, RBFOpt [4]⁴⁴4https://projects.coin-or.org/RBFOpt, StepDIRECT-0, and StepDIRECT. For both StepDIRECT-0 and StepDIRECT in Algorithm 1, we set the maximum number of function evaluations $m^{max}=2000$. The same function evaluation limit is used for DIRECT. In Algorithm 2, we select $\delta=1,\delta_{min}=0.001,\delta_{max}=2.5$, $\tau=1.5,|{\cal{D}}|=5$, and $t^{max}=1.5p$. For all other algorithms, we use the default settings. We run all algorithms for 20 times and report the mean and standard deviation results for final objective function values (denoted by “obj”) and running times in seconds (denoted by “time”) in Table 2.

From Table 2, we see that the function values returned by StepDIRECT-0 are better than those of the original DIRECT. It illustrates the benefit of our proposed strategy for identifying potentially optimal hyper-rectangles and dividing hyper-rectangles which efficiently exploits the stepwise function structure. For DIRECT and StepDIRECT-0, their objective function outputs do not change for different runs since they are deterministic algorithms.

To further see the advantage of local search incorporated into StepDIRECT-0, in StepDIRECT the local search is initialized with the best point in the potentially optimal hyper-rectangle and runs with search directions $\mathcal{D}$ randomly generated by the coordinate strategy. Compared with StepDIRECT-0, we notice that the StepDIRECT algorithm achieved lower objective function values. From Table 2, StepDIRECT shows the best overall performance in terms of solution quality except for mpg dataset. By embedding the local search, we can significantly improve the solution quality. In general, the proposed algorithm runs faster than other baseline methods DE, PSO, and RBFOpt, except for the run time of bodyfat dataset, DE outperforms StepDIRECT.

We tune the hyper-parameters for: 1) multi-class classification with linear support vector machines (SVM) and logistic regression, and 2) imbalanced binary classification with RBF-kernel SVM [11]. We use three datasets: MNIST from [18], PenDigits from the UCI Machine Learning Repository [6], and a synthetic dataset Synth generated by the Mldatagen generator [29]. For all datasets, we set the ratio among training, validation and test data partitions as $3:1:1$ and report the performance on the test dataset.