Stochastic optimization with adaptive restart: a framework for integrated local and global learning

Abstract

A common approach to global optimization is to combine local optimization methods with random restarts. Restarts have been used as a performance boosting approach. They can be a means to avoid “slow progress” by exploiting a potentially good solution, and restarts can enable the potential discovery of multiple local solutions, thus improving the overall quality of the returned solution. A multi-start method is a way to integrate local and global approaches; where the global search itself can be used to restart a local search. Bayesian optimization methods aim to find global optima of functions that can only be point-wise evaluated by means of a possibly expensive oracle. We propose the stochastic optimization with adaptive restart (SOAR) framework, that uses the predictive capability of Gaussian process models as a means to adaptively restart local search and intelligently select restart locations with current information. This approach attempts to balance exploitation with exploration of the solution space. We study the asymptotic convergence of SOAR to a global optimum, and empirically evaluate SOAR performance through a specific implementation that uses the Trust Region method as the local search component. Numerical experiments show that the proposed algorithm outperforms existing methodologies over a suite of test problems of varying problem dimension with a finite budget of function evaluations.

Appendices

Appendix A: details of Gaussian processes

The basic idea behind meta-models such as Gaussian processes is that any smooth function $f(\textit{\textbf{x}})$ can be interpreted as a realization from a stationary Gaussian process. Given n points in $\mathbb {X}$, $\mathbb {S}=\{\textit{\textbf{x}}_1, \ldots , \textit{\textbf{x}}_n\}$, and associated function values, $f(\textit{\textbf{x}}_i$), for $i=1, \ldots , n$, we construct a stationary Gaussian process $F(\textit{\textbf{x}})=\mu + Z(\textit{\textbf{x}})$ around the smooth function $f(\textit{\textbf{x}})$ using the n observed points. Here, $\mu$ is the constant mean (more complex models can be considered  [44]), and $Z(\textit{\textbf{x}})$ is the Gaussian process $Z(\textit{\textbf{x}})\sim GP(0,\tau ^2\mathbf {R})$ with spatial correlation matrix $\mathbf {R}$ and overall process variance $\tau ^2$.

We adopt the Gaussian correlation function such that $R_{ij} = \prod _{l=1}^d e^{-(\theta _l|x_{il}-x_{jl}|)^2}$, for $i,j = 1, \ldots , n$, and where d is the dimension of the vector of hyper parameters $\varvec{\theta }$ that collects the correlation factors in charge of smoothing the predictor with varying intensity. The two parameters $\mu$ and $\tau ^2$ can be estimated through the following maximum likelihood estimators  [44],

$$\begin{aligned} \hat{\mu }&= \frac{\textit{\textbf{1}}_n^T\mathbf {R}^{-1}\textit{\textbf{f}}}{\textit{\textbf{1}}^T_n\mathbf {R}^{-1}\textit{\textbf{1}}_n}, \\ \tau ^2&= \frac{(\textit{\textbf{f}}-\textit{\textbf{1}}_n\hat{\mu })^T \mathbf {R}^{-1}(\textit{\textbf{f}}-\textit{\textbf{1}}_n\hat{\mu })}{n} \end{aligned}$$

where $\textit{\textbf{f}}$ represents the n-dimensional vector of the function evaluations at the sampled points within the set $\mathbb {S}=\left\{ \textit{\textbf{x}}_{1},\textit{\textbf{x}}_{2},\ldots ,\textit{\textbf{x}}_{n}\right\}$, and $\textit{\textbf{1}}_n$ is an n-vector of ones.

Thus, we can build a predictor function $\hat{f}(\textit{\textbf{x}})$ and model variance $\hat{s}^2(\textit{\textbf{x}})$ for $\textit{\textbf{x}}\in \mathbb {X}$ following  [44].

The predictor function is defined for any $\textit{\textbf{x}} \in {\mathbb {X}}$:

$$\begin{aligned} \hat{f}(\textit{\textbf{x}}) = \hat{\mu } + \mathbf {r}(\textit{\textbf{x}})^T\mathbf {R}^{-1}(\textit{\textbf{f}} - \mathbf {1}_n\hat{\mu }) \end{aligned}$$

with a corresponding predictive model variance of:

$$\begin{aligned} \text{ Var }[\hat{f}(\textit{\textbf{x}})]= & {} \hat{s}^2(\textit{\textbf{x}})\\= & {} \tau ^2\left( 1-\mathbf {r}(\textit{\textbf{x}})^T\mathbf {R}^{-1}\mathbf {r}(\textit{\textbf{x}}) + \frac{(1-\mathbf {1}_n^T\mathbf {R}^{-1}\mathbf {r}(\textit{\textbf{x}}))^2}{\mathbf {1}_n^T\mathbf {R}^{-1}\mathbf {1}_n} \right) . \end{aligned}$$

where $\mathbf {r}$ is the n-dimensional vector of the correlations between the predicted variance at location $\textit{\textbf{x}}$ and the observed model error at the sampled locations $\textit{\textbf{x}}\in \mathbb {S}$:

$$\begin{aligned} r_i(\textit{\textbf{x}}) = \text{ Corr }(Z(\textit{\textbf{x}}),Z(\textit{\textbf{x}}_i)),\textit{\textbf{x}}_i\in \mathbb {S}. \end{aligned}$$

The interested reader can refer to  [44] and references.

Appendix B: details of the trust region algorithm

The implementation of SOAR in this paper uses the trust region method  [11, 33] as its local search, an iterative derivative-free optimization approach that, under appropriate conditions, is guaranteed to converge to a local minimum  [33]. An iteration of the trust region algorithm approximates the objective function (adopting usually a linear or quadratic model) around a location referred to as centroid of the search, and proceeds to minimize the approximating function over a “trust region” (usually a simplex or ellipsoid). The approximation allows to easily solve the minimization problem. The solution is used as the centroid location for the next iteration if the linear/quadratic surface is a good approximation for the true function within the trust region.

More specifically, the $\ell _{k}^\mathrm{th}$ iteration of our implementation of the trust region method starts with a centroid $\tilde{\textit{\textbf{x}}}^{c}_{\ell _{k}}$ of an associated hypercube as a trust region, denoted $\mathbb {R}_{\ell _{k}}$, with length of the side of the hypercube equal to $r_{\ell _{k}}$. An iteration of the algorithm results in a centroid update and/or a trust region size update, achieved by minimizing a quadratic model, $\hat{f}_Q(\textit{\textbf{x}})$. The quadratic model is given by:

$$\begin{aligned} \hat{f}_Q({\textit{\textbf{x}}}) =&f(\varvec{\tilde{x}}^{c}_{\ell _{k}}) + \nabla f(\varvec{\tilde{x}}^{c}_{\ell _{k}})(\varvec{\tilde{x}}^{c}_{\ell _{k}} - \textit{\textbf{x}}) \nonumber \\&+ \frac{1}{2}\nabla ^2 f(\varvec{\tilde{x}}^{c}_{\ell _{k}})(\varvec{\tilde{x}}^{c}_{\ell _{k}} - \textit{\textbf{x}})^2 \end{aligned}$$

(13)

where $\nabla f(\varvec{\tilde{x}}^{c}_{\ell _{k}})$ is an approximation of the gradient at the centroid, and $\nabla ^2 f(\varvec{\tilde{x}}^{c}_{\ell _{k}})$ is an approximation of the Hessian matrix. Hence, the Trust Region Sub-problem (TRS):

$$\begin{aligned} \text{(TRS) } \ {\tilde{\textit{\textbf{x}}}}^{*}_{\ell _{k}}\arg \min&\ \ \hat{f}_Q({\textit{\textbf{x}}}) \nonumber \\ s.t.\&\ \ {\textit{\textbf{x}}} \in \mathbb {R}_{\ell _{k}}. \end{aligned}$$

(14)

where ${\tilde{\textit{\textbf{x}}}}^{*}_{\ell _{k}}$ becomes the candidate centroid for the next trust region iteration.

Two decisions need to be made by the algorithm to progress to the next iteration: (1) whether to accept the candidate centroid ${\tilde{\textit{\textbf{x}}}}^{*}_{\ell _{k}}$; and (2) whether to shrink or expand the size parameter $r_{\ell _{k}}$. The ratio-comparison test answers these questions by comparing the true function value $f(\tilde{\textit{\textbf{x}}}^{*}_{\ell _{k}})$ to the predicted value from the quadratic model $\hat{f}_Q(\tilde{\textit{\textbf{x}}}^{*}_{\ell _{k}})$, and constructing the following statistic  [9, 33],

$$\begin{aligned} \rho = \frac{f(\varvec{\tilde{x}}^{c}_{\ell _{k}})-f(\tilde{\textit{\textbf{x}}}^{*}_{\ell _{k}})}{\hat{f}_Q(\varvec{\tilde{x}}^{c}_{\ell _{k}})-\hat{f}_Q(\tilde{\textit{\textbf{x}}}^{*}_{\ell _{k}})}. \end{aligned}$$

(15)

when $\rho > 1$, then we have produced a better than predicted reduction over the objective function, while $\rho < 0$ indicates that the current model is inadequate over the current trust region, and while a reduction in the objective value was predicted, the trust region step produced a worsening solution, i.e., $f(\varvec{\tilde{x}}^{c}_{\ell _{k}})<f(\tilde{\textit{\textbf{x}}}^{*}_{\ell _{k}})$.

Comparing $\rho$ with the user defined threshold values $\eta _{1},\eta _{2}$, $0<\eta _1\le \eta _2<1$, we have three possible cases:

Case 1:

$\rho \le \eta _1 \implies$ keep current centroid $\tilde{\textit{\textbf{x}}}^{c}_{\ell _{k}+1} \leftarrow \tilde{\textit{\textbf{x}}}^{c}_{\ell _{k}}$; reduce trust region size $r_{\ell _{k}+1} \leftarrow r_{k_\ell } \cdot \omega$.

Case 2:

$\eta _1 < \rho \le \eta _2 \implies$ accept candidate centroid $\tilde{\textit{\textbf{x}}}^{c}_{\ell _{k}+1} \leftarrow \varvec{\tilde{x}}^{*}_{\ell _{k}}$; keep current trust region size $r_{\ell _{k}+1}\leftarrow r_{\ell _{k}}$.

Case 3:

$\rho > \eta _2 \implies$ accept candidate centroid $\tilde{\textit{\textbf{x}}}^{c}_{\ell _{k}+1} \leftarrow \varvec{\tilde{x}}^{*}_{\ell _{k}}$; expand trust region size $r_{\ell _{k}+1} \leftarrow r_{\ell _{k}} \cdot \gamma$.

where $\omega \in (0,1)$ is the trust region reduction rate and $\gamma > 1$ is the trust region expansion rate.

Under Case 1, we reformulate the subproblem (TRS) with the same local quadratic model but under more restrictive trust region constraints. Under Case 2 or Case 3, we recompute $\nabla f(\tilde{\textit{\textbf{x}}}^{c}_{\ell _{k}})$ and $\nabla ^2 f(\tilde{\textit{\textbf{x}}}^{c}_{\ell _{k}})$ in order to update the quadratic model $\hat{f}_Q(\textit{\textbf{x}})$, since the trust region centroid has moved.

The following stopping criteria are used (in addition to the restart-related conditions):

Criterion 1:

$\Vert \nabla f(\tilde{\textit{\textbf{x}}}^{*}_{\ell _{k}}) \Vert < \epsilon _1$

Criterion 2:

$\Vert \varvec{\delta }_{\ell _{k}} \Vert < \epsilon _2$

where $\Vert \cdot \Vert$ is the Euclidean norm, $\epsilon _{1}, \epsilon _{2} > 0$ are user-defined parameters, and $\varvec{\delta }_{\ell _{k}} = \tilde{\textit{\textbf{x}}}^{*}_{\ell _{k}} - \tilde{\textit{\textbf{x}}}^{c}_{\ell _{k}}$ represents the step size between the current and proposed centroid. Note that $\mathbb {R}_{\ell _{k}} \rightarrow \mathbb {X}$ as $\ell _{k} \rightarrow \infty$, and $||\varvec{\delta }_{\ell _{k}}||_{2} \rightarrow 0$ as $\ell _{k} \rightarrow \infty$.

It was shown in  [33] (Theorem 4.8) that, when the iterations grow $\ell _{k}\rightarrow \infty$:

$$\begin{aligned} \tilde{\textit{\textbf{x}}}^{c}_{\infty } \in {\mathbb {X}} : f(\tilde{\textit{\textbf{x}}}^{c}_{\infty }) \le f(\textit{\textbf{x}}), \ \forall \textit{\textbf{x}} \in \mathcal {N}_\epsilon (\tilde{\textit{\textbf{x}}}^{c}_{\infty }), \ \forall \ \epsilon > 0 \end{aligned}$$

(16)

regardless of the initial starting location of the trust region search, where $\tilde{\textit{\textbf{x}}}^{c}_{\infty }$ is the centroid in the limit, and $\mathcal {N}_\epsilon (\tilde{\textit{\textbf{x}}}^{c}_{\infty })$ is an $\epsilon$ radius ball centered at $\tilde{\textit{\textbf{x}}}^{c}_{\infty }$. The result requires (1) Lipschitz continuity for f, (2) boundedness of f, (3) $\eta _1 \in (0,0.25)$, and (4) uniform boundedness in norm of the Hessian of f by a constant $\beta$.