A derivative-free optimization algorithm for the efficient minimization of functions obtained via statistical averaging

Abstract

This paper considers the efficient minimization of the infinite time average of a stationary ergodic process in the space of a handful of design parameters which affect it. Problems of this class, derived from physical or numerical experiments which are sometimes expensive to perform, are ubiquitous in engineering applications. In such problems, any given function evaluation, determined with finite sampling, is associated with a quantifiable amount of uncertainty, which may be reduced via additional sampling. The present paper proposes a new optimization algorithm to adjust the amount of sampling associated with each function evaluation, making function evaluations more accurate (and, thus, more expensive), as required, as convergence is approached. The work builds on our algorithm for Delaunay-based Derivative-free Optimization via Global Surrogates (${\varDelta }$-DOGS, see JOGO https://doi.org/10.1007/s10898-015-0384-2). The new algorithm, dubbed $\alpha$-DOGS, substantially reduces the overall cost of the optimization process for problems of this important class. Further, under certain well-defined conditions, rigorous proof of convergence to the global minimum of the problem considered is established.

Appendices

Appendix 1: Polyharmonic spline regression

The algorithm described in this paper depends upon a smooth regression $p^k(x)$ (see Assumption 1). The best technique for computing the regression is problem dependent. As with [12,13,14], a key advantage of our Delaunay-based approach in the present work is that it facilitates the use of any suitable regression technique, subject to it satisfying the “strict” regression property given in Definition 4. Since our numerical tests all implement the polyharmonic spline regression technique, the derivation of this regression technique is briefly explained in this appendix; additional details may be found in [46].

The polyharmonic spline regression p(x) of a function f(x) in ${\mathbb {R}}^n$ is defined as a weighted sum of a set of radial basis functions $\varphi (r)$ built around the location of each measurement point, plus a linear function of x:

$$\begin{aligned}&p(x) = \sum _{i = 1}^N w_i\,\varphi (r) + v^T \begin{bmatrix} 1 \\ x \end{bmatrix}, \nonumber \\&\quad \text {where} \quad \varphi (r) = r^3 \quad \text {and} \quad r = {\Vert }x - x_i{\Vert }. \end{aligned}$$

(47)

The weights $w_i$ and $v_i$ represent N and $n+1$ unknowns. Assume that $\{y(x_1), y(x_2), \dots , y(x_n)\}$ is the set of measurements, with standard deviations $\{\sigma _1, \sigma _2, \dots ,\sigma _2 \}$. The $w_i$ and $v_i$ coefficients are computed by minimizing the following objective function, which expresses is a tradeoff between the fit to the observed data and the smoothness of the regressor:

$$\begin{aligned} L_p(x)= \sum _{i=1}^N \left[ \frac{(p(x_i)-y(x_i))}{\sigma _i} \right] ^2 +\lambda \int _{B} {{|}\nabla ^m p(x){|}}, \end{aligned}$$

(48)

where B is a large box domain containing all of the $x_i$, and $\nabla ^m p(x)$ is the vector including all m derivatives of p(x) (see [20]). It is shown in [20] that the first-order optimality condition for the objective function (48) is as follows:

$$\begin{aligned} p(x_i)-y(x_i)+ \rho \,\sigma _i^2 w_i=0, \quad \forall 1 \le i \le N, \end{aligned}$$

(49)

where $\rho$ is a parameter proportional to $\lambda$. In summary, the coefficient of the regression can be derived by solving:

$$\begin{aligned}&\begin{bmatrix} F &{} V^T \\ V &{} 0 \end{bmatrix} \begin{bmatrix} w \\ v \end{bmatrix}= \begin{bmatrix} f(x_i) \\ 0 \end{bmatrix}, \nonumber \\&\quad F_{ij} = \varphi ({\Vert }x_i-x_j{\Vert }) +\rho \delta _{i,j} \,\sigma _i^2, \quad \quad V = \begin{bmatrix} 1 &{} 1 &{} \dots &{} 1 \\ x_1 &{} x_2 &{} \dots &{} x_N \end{bmatrix}, \end{aligned}$$

(50)

where $\delta _{i,j}$ is the Kronecker delta.

The problem which is left to solve when computing the regression is to find an appropriate value of $\rho \in [0, \infty )$. Solving (50) for any value of $\rho$ gives a unique regression, denoted $p(x,\rho )$. The parameter $\rho$ is then obtained by a predictive mean-square error criteria developed in §4.4 in [46], which is given by imposing the following condition:

$$\begin{aligned} T(\rho )=\sum _{i=1}^N \left[ \frac{p(x_i,\rho )-y(x_i)}{\sigma _i}\right] ^2=1. \end{aligned}$$

(51)

For $\rho \rightarrow \infty$, $w_i\rightarrow 0$, and the solution of (50) is a weighted mean-square linear regression, which is obtained by solving (51). If $T(\infty ) \le 1$, we take this linear regression as the best current regression for the available data. Otherwise, we have $T(\infty )>1$ and (by construction) $T(0)=0$; thus, (51) has a solution with finite $\rho >0$, which gives the desired regression.

If $T(\infty ) > 1$, we thus seek a $\rho$ for which for $T(\rho )=1$. Following [46], using (50), (51) simplifies to:

$$\begin{aligned} T(\rho )=\rho ^2 \left( \sum _{i=1}^N w_{i}) \,\sigma _i\right) ^2 =1, \end{aligned}$$

(52)

where $w_{i, \rho }$ is the $w_i$ which is obtained by solving (50). Define Dw and Dv as the vectors whose i-th elements are the derivatives of $w_i$ and $v_i$ with respect to $\rho$, then

$$\begin{aligned}&T'(\rho )= \rho ^2 \sum _{i=1}^N w_i \, Dw_i \sigma _i^2+ 2\, \rho \left( \sum _{i=1}^N w_{i, \rho } \sigma _i\right) ^2, \\&\begin{bmatrix} F &{} V^T \\ V &{} 0 \end{bmatrix} \begin{bmatrix} Dw \\ Dv \end{bmatrix}+ \begin{bmatrix} \rho {\varSigma }_2 &{} 0 \\ 0 &{} 0 \end{bmatrix} \begin{bmatrix} w \\ v \end{bmatrix}= \begin{bmatrix} 0 \\ 0 \end{bmatrix}, \end{aligned}$$

where ${\varSigma }_2$ is a diagonal matrix whose i-the diagonal element is $\rho \,\sigma _i^2$. Therefore, the analytic expression for the derivative of $T(\rho )$ is available. Thus, (51) can be solved quickly using Newton’s method.

The regression process presented here, imposing (1) as suggested by [46], is designed to obtain a regression which is reasonably smooth. However, there is no guarantee that this particular regression satisfies the strictness property required in the present work (see Definition 4). Note, however, that by imposing $\rho =0$, the regression is made strict for arbitrary small $\beta$. Thus, to satisfy strictness for a given finite $\beta$, the value of $\rho$ must sometimes be decreased from that which satisfies (1), as necessary.

Appendix 2: UQ for finite-time-averaging of the Lorenz system

This appendix summarizes briefly the simple empirical approach that is used in this paper to quantify the uncertainty of the cost function (44) when it is estimated using a finite time average.

In the method used, we simply simulated the Lorenz system (42) 30 independent times with various initial values for (X, Y, Z). The cost function was then approximated using these different simulation lengths, and the standard deviation of the estimations obtained using various simulation lengths was calculated. The simple model given by ${A}/{\sqrt{T}}$ for the uncertainty was found to fit this empirical calculation quite well, as shown in Fig. 10.

Fig. 10

Uncertainty quantification (UQ) for finite-time-average approximations of the inifinite-time-average statistic of interest in the cost function related to the Lorenz model. The uncertainty quantification model of ${A}/{\sqrt{T}}$ is found to give a very good empirical fit

Appendix 3: Application of $\alpha$-DOGS on higher dimensional problems

We now consider the challenge of increasing the dimension of the problem considered. This is rather difficult, as the algorithm developed herein (specifically, the computational cost of computing the Delaunay triangulation) scales rather poorly as the dimension of the problem is increased. In this appendix, the convergence behavior on a 5-dimensional problem is thus investigated.

The test problem that is considered here is the stochastically-obscured scaled Schwefel function (41) for $n=5$. As in Sect. 6, performance of the $\alpha$-DOGS, ${\varDelta }$-DOGS, and SMF algorithms are characterized. In this investigation, for all datapoint in ${\varDelta }$-DOGs and SMF algorithm, a fixed averaging time of $T=100$ was used. Figure 11 shows the convergence history of the regret function for each method. It is observed that the required total averaging time for the optimization problem with $\alpha$-DOGS is significantly better than the other methods considered.

To scale to even higher dimensional problems, the practical (though, perhaps, not provably globally convergent) idea of using only approximate Delaunay triangulations has been proposed, and will be considered in future work.

Fig. 11

Convergence history of optimization algorithms for (41) for $n=5$