Safe global optimization of expensive noisy black-box functions in the $\delta$-Lipschitz framework

TL;DR

This paper introduces a novel framework and algorithms for safely maximizing noisy black-box functions under Lipschitz conditions, ensuring evaluations occur only at safe points and providing theoretical guarantees and numerical validation.

Contribution

It proposes a $ ext{delta}$-Lipschitz framework and two algorithms for safe global optimization of noisy functions, addressing safety constraints and noise challenges.

Findings

Algorithms reliably identify safe regions and maximize within them.

Theoretical convergence guarantees are established for the proposed methods.

Numerical experiments confirm the effectiveness of the algorithms.

Abstract

In this paper, the problem of safe global maximization (it should not be confused with robust optimization) of expensive noisy black-box functions satisfying the Lipschitz condition is considered. The notion "safe" means that the objective function $f(x)$ during optimization should not violate a "safety" threshold, for instance, a certain a priori given value $h$ in a maximization problem. Thus, any new function evaluation (possibly corrupted by noise) must be performed at "safe points" only, namely, at points $y$ for which it is known that the objective function $f(y) > h$. The main difficulty here consists in the fact that the used optimization algorithm should ensure that the safety constraint will be satisfied at a point $y$ before evaluation of $f(y)$ will be executed. Thus, it is required both to determine the safe region $\Omega$ within the search domain~$D$ and to find the global maximum within $\Omega$. An additional difficulty consists in the fact that these problems should be solved in the presence of the noise. This paper starts with a theoretical study of the problem and it is shown that even though the objective function $f(x)$ satisfies the Lipschitz condition, traditional Lipschitz minorants and majorants cannot be used due to the presence of the noise. Then, a $\delta$-Lipschitz framework and two algorithms using it are proposed to solve the safe global maximization problem. The first method determines the safe area within the search domain and the second one executes the global maximization over the found safe region. For both methods a number of theoretical results related to their functioning and convergence is established. Finally, numerical experiments confirming the reliability of the proposed procedures are performed.

Tables

Table 1. Table 1: Test Case 1

Problem	Lipschitz Function	Interval	L	Threshold	Ref.
1	$f\left(x\right)=-\frac{1}{6}{x}^6+\frac{52}{25}{x}^5-\frac{39}{80}{x}^4-\frac{71}{10}{x}^3+\frac{79}{20}{x}^2+x-\frac{1}{10}$	[-1.5,11]	13870	2974.180	[44]
2	$f\left(x\right)=-\sin x^3-\cos x^3$	[0,6.28]	2.2	-0.800	[30]
3	$f\left(x\right)=x-\sin 3x+1$	[0,6.5]	4	1.202	[30]
4	$f\left(x\right)=\frac{x^2-5x+6}{x^2+1}$	[-5,5]	6.5	0.671	[44]

Table 2. Table 3: Test Case 3

Problem	Lipschitz Function	Interval	L	Threshold	Ref.
11	$f\left(x\right)={\sum }_{i=1}^{5}i\sin \left[\left(i+1\right)x+i\right]+3$	[-10,10]	67	-24.335	[30]
12	$f\left(x\right)=\cos x-\sin 5x+1$	[0,7]	5.951	-0.545	[30]
13	$f\left(x\right)=\{\begin{array}{cc}\cos 5x,x\le \frac{3}{2}\pi & \\ \cos x,\mathrm{otherwise}& \end{array}$	[0,18]	4.999	-0.800	[30]
14	$f\left(x\right)=\{\begin{array}{cc}\sin x,x\le \pi & \\ \sin 5x,\mathrm{otherwise}& \end{array}$	[-10,10]	4.999	-0.800	[30]
15	$f\left(x\right)=-{\sum }_{i=1}^5\cos \left[\left(i+1\right)x\right]$	[-10,10]	18.119	-4.229	[30]
16	$f\left(x\right)=x\left|\sin x\right|+6$	[-10,10]	9.632	-0.332	[30]
17	$f\left(x\right)=\left|x\sin x\right|-1.5$	[-10,10]	9.632	-0.709	[30]
18	$f\left(x\right)=\{\begin{array}{cc}\sin x,\sin x>\cos x& \\ \cos x,\mathrm{otherwise}& \end{array}$	[-10,10]	1	-0.519	[30]

Table 3. Table 4: The number of points at which the objective function has been evaluated and the number of these evaluations executed by Algorithms 1 and 2 on the three test cases

Problem	Safe Expansion		Global Maximization		Total
	Points	Evaluations	Points	Evaluations	Points	Evaluations
1	15	41	20	31	35	72
2	10	43	82	178	92	221
3	16	44	10	19	26	63
4	40	95	39	99	79	194
5	52	191	29	33	81	224
6	59	92	27	30	86	122
7	150	209	128	209	278	418
8	53	245	68	169	121	414
9	380	1002	393	403	773	1495
10	71	95	35	44	106	139
11	210	485	170	196	380	681
12	57	184	57	97	114	281
13	99	441	55	74	154	515
14	101	441	35	61	136	502
15	93	255	89	108	182	363
16	54	136	46	81	100	217
17	59	249	74	169	133	418
18	27	191	29	35	56	226

Equations

$$|f(x_1)-f(x_2)|\le L|x_1-x_2|,0<L<\infty ,x_1,x_2\in D,$$

$$x^{\ast }=\underset{x\in \Omega \subseteq D}{\mathrm{argmax}}g(x),g(x)=f(x)+\xi (x),$$

$$|\xi (x)|\le \delta ,\delta >0,$$

$$\Omega =\{x:x\in D,g(x)\ge h\},\Omega ={\cup }_{j=1}^m{\Omega }_j,{\Omega }_i\cap {\Omega }_j=\varnothing ,i\ne j.$$

$$G(x)=\{y:y=f(x)+\xi (x),\xi (x)\in [-\delta ,\delta ]\}.$$

$$\psi (x,x_1)=g(x_1)-L|x_1-x|,x\in D,$$

$$|s(x_1)-s(x_2)|\le L|x_1-x_2|+\delta ,0<L<\infty ,0<\delta <\infty ,x_1,x_2\in D.$$

$$|g(x_1)-g(x_2)|=|f(x_1)+\xi (x_1)-f(x_2)-\xi (x_2)|\le$$

$$|f(x_1)-f(x_2)|+|\xi (x_1)|+|\xi (x_2)|\le L|x_1-x_2|+2\delta .$$

$$|g(x_1)-g(x_2)|\le L|x_1-x_2|,0<L<\infty ,x_1,x_2\in D,$$

$$\phi (x,x_i)=g(x_i)-L|x_i-x|-2\delta ,x\in D,$$

$$|\phi (x^{\prime },x_i)-\phi (x^{\prime \prime },x_i)|\le L|x^{\prime }-x^{\prime \prime }|,x^{\prime },x^{\prime \prime }\in D.$$

$$|\phi (x^{\prime },x_i)-\phi (x^{\prime \prime },x_i)|=|g(x_i)-L|x_i-x^{\prime }|-2\delta -g(x_i)+L|x_i-x^{\prime \prime }|+2\delta |=$$

$$|(L|x_i-x^{\prime \prime }|-L|x_i-x^{\prime }|)|=L|(|x_i-x^{\prime \prime }-x^{\prime }+x^{\prime }|-|x_i-x^{\prime }|)|\le$$

$$L|(|x^{\prime }-x^{\prime \prime }|+|x_i-x^{\prime }|-|x_i-x^{\prime }|)|=L|x^{\prime }-x^{\prime \prime }|.$$

$${\Phi }_k(x)=\underset{1\le i\le k}{\max }\phi (x,x_i),x\in D,$$

$${\Phi }_k(x)\le g(x)=f(x)+\xi (x),\forall \xi (x)\in [-\delta ,\delta ],x\in D.$$

$${\Phi }_k(x)\le \underset{\xi (x)\in [-\delta ,\delta ]}{\min }g(x)=f(x)-\delta .$$

$${\Phi }_k(x)=\phi (x,x_j)=g(x_j)-L|x_j-x|-2\delta .$$

$$g(x_j)-g(x)\le |g(x_j)-g(x)|\le L|x_j-x|+2\delta .$$

$$g(x_j)-(f(x)-\delta )\le L|x_j-x|+2\delta .$$

$${\Phi }_k(x)=g(x_j)-L|x_j-x|-2\delta +(f(x)-\delta )-(f(x)-\delta )\le$$

$$f(x)-\delta -L|x_j-x|-2\delta +L|x_j-x|+2\delta =f(x)-\delta .$$

$${\check{x}}_k=\min \{x_i:x_i\in \Omega ,1\le i\le k\},{\hat{x}}_k=\max \{x_i:x_i\in \Omega ,1\le i\le k\},$$

$$g({\check{x}}_k)-L({\check{x}}_k-x)-2\delta =h,x<{\check{x}}_k,$$

$$g({\hat{x}}_k)-L(x-{\hat{x}}_k)-2\delta =h,x>{\hat{x}}_k,$$

$${\check{x}}_{k+1}={\check{x}}_k-L^{-1}(g({\check{x}}_k)-2\delta -h),{\hat{x}}_{k+2}={\hat{x}}_k+L^{-1}(g({\hat{x}}_k)-2\delta -h).$$

$$h=\phi ({\check{x}}_{k+1},{\check{x}}_k)\le g({\check{x}}_{k+1}),h=\phi ({\hat{x}}_{k+2},{\hat{x}}_k)\le g({\hat{x}}_{k+2}).$$

$$y_1-y_2\ge 2\delta -\sigma$$

$$x^{\ast }=\underset{x\in \hat{\Omega }\subseteq \Omega \subseteq D}{\mathrm{argmax}}g(x),g(x)=f(x)+\xi (x),$$

Full text

Keywords: Safe global optimization; expensive black-box functions; noise; Lipschitz condition; Machine learning.

1 Introduction

Many complex industrial applications are characterized by black-box, multiextremal, and expensive objective functions that should be optimized (see [1, 11, 21, 26, 31, 37, 40, 43, 46]). The word expensive means here that each evaluation of the objective function $f(x)$ is a time consuming operation. Since locally optimal solutions often do not give a sufficiently high level of improvement w.r.t. a currently available solution obtained by engineers using practical reasons, global optimization problems are considered (see [10, 1, 13, 16, 20, 21, 33, 40, 47, 51, 49, 50]). Unfortunately, a practical global optimization process is often performed under a limited budget, i.e., the number of allowed evaluations of $f(x)$ is fixed a priori and is not very high (see a detailed discussion in [43]) requiring so an accurate development of fast global optimization methods (see, e. g., [3, 2, 18, 27, 31, 32, 40, 47, 50]).

Recently, a class of real applications having an additional important constraint on the value of the objective function is under an accurate study (see, e. g., [5, 14, 48] in the context of Lipschitz optimization, control problems, reinforcement learning, Bayesian optimization, etc.). This constraint requires that the objective function $f(x)$ during optimization should not violate a “safety” threshold (that should not be confounded with robust optimization, see, e. g., [4]), for instance, a certain a priori given value $h$ in a maximization problem. Thus, any new function evaluation must be performed at “safe points” only, namely, at points $y$ for which it is known that the objective function value will not violate the safety threshold $h$, i.e., it should be $f(y)>h$. The main difficulty here consists in the fact that the used optimization algorithm should ensure that the safety constraint will be satisfied at a point $y$ before evaluation of $f(y)$ will be executed.

This requirement is very important and represents a key difference w.r.t. traditional constraint problems where constraints can be checked and the algorithm then behaves in dependence on whether a constraint was satisfied or not. In this kind of problems, it is often allowed to evaluate the objective function $f(x)$ not only at admissible but also at inadmissible points. In contrast, in the safe optimization problems evaluation of $f(x)$ at inadmissible points is strictly prohibited. Thus, while formulations using global optimization with unknown constraints is well suited for simulation-optimization problems, the safe global optimization formulation is more appropriate for online control and online learning and optimization problems (see, e. g., [5, 14, 48]).

For instance, let us consider, as an example of the safe optimization problem, tuning parameters of a controller used for automatic process control in the manufacturing industry, especially in the pneumatic, electronic, robotic and automotive domains (see, e. g., [12, 35]). One of the most widely used types of controllers is known as PID, according to the three components: Proportional, Integral, and Derivative. A PID controller implements a control loop mechanism based on feedback: depending on the error value, computed continuously as the difference between a desired setpoint and a measured process variable, the PID controller applies a correction based on its proportional, integral, and derivative terms. Although a software model of the system to control can be used to initially design the PID controller, an expensive and prone-to-failures manual tuning phase is needed to set up, on the real-life system, the optimal values of its parameters from where the need to perform this phase “safely” arises. Depending on the characteristics of the system to control, the safety of the controller is evaluated with respect to one or more of the following issues: responsiveness of the controller, overshooting of the desired setpoint by the measured process variable as well as oscillatory behaviour around the setpoint. In particular, overshooting and oscillatory behaviour could lead to malfunctioning or breakage of the system to control. Thus, searching for the optimal tuning of the PID parameters on a real system, without considering “safety” of the system itself, could damage it, in some cases irreparably.

In some sense, it is required both to determine the safe region $\Omega$ within the search domain $D$ and to find the global optimum within the safe region $\Omega$. Usually, at least one safe point is known before the start of the optimization process (e. g., it is taken using the current working configuration of the optimized industrial system). Thus, it is required to invent a “safe expansion” mechanism to extend the current safe region from the starting point in order to find the whole $\Omega$. It can happen that, if the safe region $\Omega$ consists of several disjoint subregions, those subregions which do not contain initially provided safety points will be never found. In this case, it is not possible to talk about the global optimum over the whole safe region $\Omega$ and a global optimum over the current safe subregion should be found.

A further complication that is frequently present in applied optimization problems (see [1, 13, 21, 29, 33, 47, 49, 50]) is the presence of noise affecting evaluations of the objective function $f(x)$. These evaluations should remain safe even when they are corrupted by noise. It is difficult to underestimate the importance of taking into consideration the presence of noise since it can make invalid many assumptions (e. g., convexity, derivability, Lipschitz continuity, etc.) usually done w.r.t. optimized functions. In spite of its crucial impact, noise is often not considered in detail in safe global optimization problems whereas the afore mentioned approaches [5, 14, 48] do it and provide some probability-based considerations on safety. Precisely these papers have stimulated us to study safe global optimization problems with noise.

Since Lipschitz continuity is a quite natural assumption for applied problems (specifically for technical systems, see, e. g., [17, 20, 25, 33, 40, 47]), we consider here objective functions that satisfy the Lipschitz condition over the search domain $D$. Thus, our problem becomes Lipschitz global optimization problem broadly studied in the literature (see, e. g., [3, 15, 17, 18, 21, 25, 33, 39, 40, 47] and references given therein). In this paper, we propose a new Lipschitz-based safe global optimization algorithm, specifically designed to work in the noisy setting. It is proved that, in spite of the presence of the noise, our approach does not permit any violation of the safety threshold. The only assumption made with respect to the noise is its boundedness, with a maximal level of the noise known a priori.

The remaining part of the manuscript is structured as follows. Section 2 contains statement of the problem and its analysis. Section 3 presents a theoretical investigation of a reliable expansion of the safe region and proposes an algorithm realizing this expansion. Section 4 introduces a global maximization algorithm working in the presence of the noise over the found safe region. Section 5 proposes three series of numerical experiments confirming theoretical results and showing a reliable performance of the two introduced methods. Section 6 concludes the paper.

2 Statement of the problem and its analysis

In order to start, let us present a general formulation of the safe global optimization problem using the Lipschitz framework in one dimension. As was already mentioned, Lipschitz global optimization problems can be very often encountered in applications even when $f(x)$ is univariate. Nowadays problems of this kind with and without noise are under an intensive study (see, e. g., [7, 8, 6, 9, 11, 22, 23, 24, 26, 28, 29, 33, 36, 37, 38, 45, 40, 41, 42, 44]).

To state the problem formally, let us suppose that a function $f(x)$ satisfies over a search domain $D=[a,b]$ the Lipschitz condition

[TABLE]

with an a priori known Lipschitz constant $L$111Clearly, there can exist several constants $\tilde{L}$ such that if they are placed in (1) instead of $L$ the inequality will hold. Without loss of generality we suppose hereinafter that $L>\min\tilde{L}$.. Then, given a safety threshold $h>0$, it is required, in the presence of noise ${\xi(x)}$, to find an approximation of the point $x^{*}$ and an estimate of the corresponding value $g(x^{*})$ such that

[TABLE]

where $\Omega$ is the safe region that can consist of several disjoint subregions $\Omega_{j},1\leq j\leq m,$ and the noise ${\xi(x)}$ is bounded by a known value $\delta$, i.e.,

[TABLE]

[TABLE]

Thus, at each point $x$ the value $g(x)$ can belong to the set

[TABLE]

It should be stressed that in traditional noisy optimization problems (see, e. g., [8]) the maximizer of $f(x)$, not of $g(x)$, is of interest. This is not the case in safe optimization since the function $g(x)$ is measured and not $f(x)$ and, therefore, it is important that precisely $g(x)$ would be as larger than the threshold $h$ as possible.

In order to illustrate the problem (1)–(5), let us consider Fig. 1. It shows the Lipschitz function $f(x)=sin(x)\cdot cos(4x)$ having $L=4$ and defined over $D=[-2,2]$. The safety threshold is $h=-0.8$ and the initial safe point is $x_{1}=-0.8$. It can be seen that $f(x_{1})\neq g(x_{1})\in G(x_{1})$, in particular $g(x_{1})<f(x_{1})$. The unsafe region is shown in red; it consists of two subregions. Notice that in the right-hand unsafe subregion $f(x)>h$ but due to the presence of the noise this subregion becomes unsafe. The safe region, $Q$, consists of three subregions $Q_{1},Q_{2},$ and $Q_{3}$. The presence of the safe point $x_{1}$ in the central subregion $Q_{2}$ and the knowledge of the Lipschitz constant $L$ together with supposition $g(x_{1})<f(x_{1})$ allow us to build the following Lipschitz minorant

[TABLE]

and to determine the initial safety region shown in green. It should be expanded during optimization using the Lipschitz information, however, since we have only one initial safe point at $Q_{2}$, only this safe subregion can be explored and the left-hand and right-hand safe subdomains $Q_{1}$ and $Q_{3}$ will remain undiscovered.

A strategy for expanding the initial safe region is illustrated in Fig. 2. Clearly, the maximal possible expansion of the currently found safe region can be obtained by evaluating the objective function at the points $x_{2}$ and $x_{3}$ being its extrema. Unfortunately, if the noise is not taken into account appropriately, there is a risk to overestimate the safe region. In real-life applications, at each safe point $y$ we do not know whether $g(y)\leq f(y)$ or $g(y)>f(y)$. For instance, it can be seen in Fig. 2 that in this example both $g(x_{2})>f(x_{2})$ and $g(x_{3})>f(x_{3})$. Therefore, functions $\psi(x,x_{2})$ and $\psi(x,x_{3})$ constructed using the Lipschitz constant $L$ of $f(x)$ at the points $x_{2}$ and $x_{3}$ are not minorants for $g(x)$ anymore. In particular, this fact results in an error in determining the safe region starting from the point $x_{2}$. It can be seen that the region shown in Fig. 2 in hatched red is unsafe but it is considered by the described procedure of the expansion to be safe. The source of this error is the fact that even though the objective function $f(x)$ satisfies the Lipschitz condition (1), the function $g(x)$ is not Lipschitzian as it is shown in Theorem 1 below. To prove it, we need the following definition (notice that similar but slightly different functions have been considered in [49]).

Definition 1**.**

A function $s(x)$ is called $\delta$-Lipschitz over the interval $D$ if it satisfies the following condition

[TABLE]

Theorem 1**.**

Suppose that the function $f(x)$ satisfies the Lipschitz condition (1) with a constant $L$, then the following two assertions hold:

i. The function $g(x)$ from (2) does not satisfy the Lipschitz condition.

ii. The function $g(x)$ from (2) is $2\delta$-Lipschitzian where $\delta$ is from (3).

Proof.

Let us first prove the second assertion. It follows from definition of $g(x)$ in (2), from the Lipschitz condition (1), and from the boundedness of noise in (3) that

[TABLE]

[TABLE]

Let us now consider the assertion (i). If inequality

[TABLE]

was satisfied (where $L$ is from (1)), then it would be that $\lim_{x_{1}\rightarrow x_{2}}|g(x_{1})-g(x_{2})|=0$. However, due to (7), $\lim_{x_{1}\rightarrow x_{2}}|g(x_{1})-g(x_{2})|=2\delta$ and, therefore, $g(x)$ does not satisfy (8). ∎

3 A reliable expansion of the safe

region

In order to provide a reliable expansion from the initial safe region, it is necessary to ensure that at a new point $z$ chosen to evaluate $g(x)$ condition (4) holds, i.e., $g(z)\geq h$. Moreover, this condition should be satisfied for any value of noise from (3) and, as a consequence, for any value of $g(x)\in G(x)$ from (5). Let us start to introduce such a mechanism by supposing that $g(x)$ has been evaluated at several safe points $x_{i}\in\Omega,1\leq i\leq k,$ and introduce the following functions

[TABLE]

where $\delta$ is from (3). Then, the following result holds.

Lemma 1**.**

Function $\varphi(x,x_{i})$ satisfies Lipschitz condition (1).

Proof.

To prove Lemma 1 it is required to show that

[TABLE]

For any arbitrary $x^{\prime},x^{\prime\prime}\in D$ it follows that

[TABLE]

[TABLE]

[TABLE]

Thus, Lemma has been proven. ∎

Theorem 2**.**

Let us construct the following function

[TABLE]

then for any value of noise $\xi(x)\in[-\delta,\delta]$ it follows that $\Phi_{k}(x)$ is a minorant for the function $g(x)$ satisfying condition (7) over the search region $D$, namely, it follows

[TABLE]

Proof.

Let us consider an arbitrary point $x\in D$. In order to prove (11), it is sufficient to show that

[TABLE]

It follows from (9) that there exists an index $j,\,1\leq j\leq k,$ such that

[TABLE]

Then, since $g(x)$ satisfies condition (7), we obtain

[TABLE]

Notice that this estimate holds for any value of $g(x)\in G(x)$ including, therefore, the minimal possible value of $g(x)=f(x)-\delta$. As a result we have

[TABLE]

By adding and subtracting $f(x)-\delta$ to (12) and combining it with (13) we get

[TABLE]

[TABLE]

Since the point $x$ has been chosen arbitrarily, the last inequality concludes the proof. ∎

Corollary 1**.**

Function $\Phi_{k}(x)$ from (10) satisfies Lipschitz condition (1).

Proof.

Due to Lemma 1, functions $\varphi(x,x_{i}),1\leq i\leq k,$ satisfy Lipschitz condition (1). Since $\Phi_{k}(x)$ is maximum of $k$ Lipschitz functions satisfying (1), this condition holds for $\Phi_{k}(x)$, as well. ∎

Theorems 1 and 2 allow us to substitute unsafe Lipschitz functions $\psi(x,x_{i})$ with safe $\delta$-Lipschitz minorants $\varphi(x,x_{i})$ and give a direct suggestion on how the initial safe region should be expanded. The process of expansion is illustrated in Fig. 3. We evaluate $g(x)$ at the initial safe point $x_{1}$ (the result of evaluation $g(x_{1})$ is shown in Fig. 3 by sign “x”) and construct the minorant $\varphi(x,x_{1})$ (the value $\varphi(x_{1},x_{1})=g(x_{1})-2\delta$ is shown by black dot). Then, by using $\varphi(x,x_{1})$, we obtain the initial safe subregion $[x_{4},x_{5}]$. Clearly, since $\varphi(x,x_{1})\leq\psi(x,x_{1})$, it follows (cf. Fig. 2) that $[x_{4},x_{5}]\subset[x_{2},x_{3}]$. To expand the safe region we evaluate $g(x_{4})$ and $g(x_{5})$ and construct $\varphi(x,x_{4})$ and $\varphi(x,x_{5})$. Theorem 2 ensures that the obtained region $[x_{6},x_{7}]$ (shown in green) is safe. The process is repeated and evaluations of $g(x)$ are iteratively performed at the boundaries of the currently found safe region as follows.

To continue our analysis let us suppose for simplicity that $\Omega$ is a simply connected region, i.e., $m=1$ in (4) (if this is not the case then considerations analogous to what follows should be executed for all zones $\Omega_{j},1\leq j\leq m$). Assume that $g(x)$ has been evaluated at safe points $x_{i},1\leq i\leq k,$ and let us indicate the points

[TABLE]

being the current minimal and maximal safe points, respectively. Thus, the current found safe region is $[\check{x}_{k},\hat{x}_{k}]$. By using functions $\varphi(x,\check{x}_{k})$ and $\varphi(x,\hat{x}_{k})$ new safe points $\check{x}_{k+1}$ and $\hat{x}_{k+2}$ can be obtained as solutions of equations

[TABLE]

[TABLE]

and calculated as follows

[TABLE]

It can be easily seen that the points $\check{x}_{k+1}$ and $\hat{x}_{k+2}$ are indeed safe. In fact, it follows from Theorem 2 that due to the construction

[TABLE]

Thus, the introduced $\delta$-Lipschitz framework allows us to construct a safe expansion mechanism based on the knowledge of $L$ and $\delta$ and to obtain estimates of the left and right margins of the safe region $\Omega$ (let us call them $l$ and $r$, respectively). However, there exists an additional trouble related to situations taking place due to the presence of the noise when the process of expansion approaches borders of the current safe region. Let us illustrate this difficulty by returning once again to Fig. 3. It can be seen that the value $\varphi(x_{4},x_{4})$ is very close to the lower boundary of $G(x_{4})$ equal to $f(x_{4})-\delta$ because $g(x_{4})$ is very close to the upper boundary of $G(x_{4})$, namely, to $f(x_{4})+\delta$. In contrast, $\varphi(x_{5},x_{5})$ is very far away from the lower boundary of $G(x_{5})$ equal to $f(x_{5})-\delta$ because $g(x_{5})$ is very close to the lower boundary of $G(x_{5})$, i.e., to $f(x_{5})-\delta$. Thus, the tightness of the minorant $\varphi(x,x_{i})$ depends on the level and sign of noise $\xi(x)$. If $\xi(x)\rightarrow\delta$, then the minorant is good, if $\xi(x)\rightarrow-\delta$, then the minorant is not accurate. This fact can become crucial on the borders of the current safe region. For instance, if the next function evaluation (see Fig. 3) is executed at the point $x_{6}$ and $g(x_{6})\approx f(x_{6})+\delta$ then $\varphi(x_{6},x_{6})\approx f(x_{6})-\delta$, the constraint $\varphi(x_{6},x_{6})>h$ will be satisfied an it will be possible to execute an additional (though small) safe step on the left from $x_{6}$ and to enlarge the safe region. In contrast, if $g(x_{6})\approx f(x_{6})-\delta$ then $\varphi(x_{6},x_{6})=g(x_{6})-2\delta$ will be less than $h$ and, therefore, it will not be possible to expand the current safe region.

This observation suggests that if at a current point, $z$, that is tested for expanding the safe region using the function $\varphi(x,y),y\in\Omega,$ it follows that $\varphi(z,y)<h$, then it is necessary to continue to evaluate $g(z)$ several times trying to obtain a value close to $f(z)+\delta$ and, therefore, to obtain the best possible minorant. This strategy is illustrated in Fig. 4 where the repeated evaluations of $g(x)$ at the borders of the current safe region can be clearly seen. These repetitions allow us to establish the borders with a high precision.

In practice, the number of repetitions allowing us to obtain an expansion can be very high and closer we are to the border of the safe region smaller the expansion steps will become. Thus, it is necessary to introduce a parameter $\nu$ fixing the maximal allowed number of repetitions of evaluations of $g(x)$ at a point $z$ during the search of the border and a tolerance for the expansion process. Recall that we know that the level of noise is bounded by $\delta$. Thus, if during iterated evaluations of $g(z)$ two values, $y_{1}$ and $y_{2}$, $y_{1}>y_{2},$ have been obtained as results of these evaluations then the following condition using an a priory given accuracy $\sigma>0$ for stopping iterated evaluations of $g(x)$

[TABLE]

can be used.

Let us illustrate the work of the parameters $\nu$ and $\sigma$ by the following examples. In the first of them shown in Fig. 4 the safe region for $h=-0.85$ has been identified using parameters $\nu=15$ and $\sigma=0.05\cdot 2\delta$. The search of the safe region has been stopped after 15 repetitions at the lower and 15 repetitions at the upper margin of the current safe region, i.e., the criterion on $\nu$ has stopped the search.

In Fig. 5, the same values of $h$ and $\nu$ have been used but $\sigma=0.1\cdot 2\delta$ has been taken. Function evaluations at the lower and upper boundaries of the safe region were, respectively, 13 and 10. Namely, at both sides, the search was stopped due to the criterion on $\sigma$. Notice that the function $f(x)$ is the same in both experiments shown in Figs. 4 and 5. However, due to the presence of the noise the resulting functions $g(x)$ are different.

It should be stressed that, even if the objective function evaluations are repeatedly performed at the boundaries of the estimated safe region to better define them, there is no way to go through unsafe zones to other safe subregions (that could contain the desired global maximum) if these subregions do not contain originally provided safe points. In fact, in this specific example the safe subregion obtained starting from the initial safe point $x=0.5$ shown in Fig. 5 in green does not coincide with the entire safe region $\Omega$. As a result, the identified safe subregion does not contain the global maximum over the whole $\Omega$ since it is located within the unfound subregion of $\Omega$.

The importance of the strategy consisting in insisting on repeated evaluations on the borders of the current search region is shown even better in Fig. 6 where the same function $f(x)$ is considered but the safety threshold is decreased to $h=-1.1$. In this case, the threshold $h$ limits the safe region on the left-hand side only, whereas the right border of $\Omega$ coincides with the right margin of the search space $D$.

It can be seen from Fig. 6 that the safe expansion mechanism is able to correctly estimate the part of the safe region $\Omega$ containing the global optimum. This happens in spite of the fact that the initial safe point (shown in Fig. 6 in green) is located quite far from the upper margin of $\Omega$ and a local minimum of $f(x)$, in between, is quite close to the safety threshold $h$. At a neighborhood of this local minimizer a huge number of function evaluations is performed trying to observe a function value which would guarantee a reliable expansion of the current safe region. Since in these subregion $G(x)$ is very close to $h$, the size of each new step of the expansion is quite small. In this example, the safe subregion has been identified using $\nu=15$ and $\sigma=0.1\cdot 2\delta$. The executed function evaluations at the lower boundary of the safe region were 6 (i.e., the process of the expansion on this side was stopped due to the criterion on $\sigma$). The expansion on the upper boundary of $\Omega$ was stopped because the right margin of the search region $D$ has been reached.

A detailed description of Algorithm 1 executing the described above procedure of the reliable expansion of the safe region can be found in Appendix.

4 Global maximization over the found safe region

Let us suppose now that the $\delta$-Lipschitz expansion procedure described in the previous section has stopped after evaluating $g(x)$ at $k$ trial points and a desired approximation $\widehat{\Omega}_{k}$ (indicated just $\widehat{\Omega}$ hereinafter) of the safe region $\Omega$ has been obtained. After that, the second phase of the algorithm, namely, global maximization over the found safe region $\widehat{\Omega}$, starts. For simplicity we suppose here, as was done in the previous section, that $\widehat{\Omega}$ is a simply connected region (all the subsequent considerations in this section can easily be extended to the case of disjoint subregions composing $\widehat{\Omega}$). Notice also that the maximization problem (2) over $\Omega$ is substituted by its approximation, i.e., by the problem

[TABLE]

The resulting problem (1), (17), (3)–(5) is solved by using all previously executed evaluations of $g(x)$ at the points $x_{i}\in\widehat{\Omega},1\leq i\leq k$ (see Algorithm 1 in Appendix). The maximization procedure constructs a majorant $\Gamma_{k}(x),x\in\widehat{\Omega},$ for $g(x)$ similarly to the construction of the minorant $\Phi_{k}(x)$ from (10). Recall that during the safe expansion procedure there could be produced points $x_{i}$ where $g(x)$ has been evaluated $\nu(x_{i})$ times (where $\nu(x_{i})$ in any case is not superior to the maximal allowed number, $\nu$, of repeated evaluations of $g(x)$ at a single point). This means that different values $g_{j}(x_{i})$ could be obtained during the $j$th evaluation, $1\leq j\leq\nu(x_{i})\leq\nu$ of $g(x)$ at a point $x_{i}$, $1\leq i\leq k$. Recall, that it was necessary to use values

[TABLE]

in order to build a better minorant $\Phi_{k}(x)$. Analogously, in order to build a better majorant $\Gamma_{k}(x)$, we shall use values

[TABLE]

Then, the function $\Gamma_{k}(x)$ is built as follows

[TABLE]

[TABLE]

where $\delta$ is from (3). Both functions, $\Phi_{k}(x)$ and $\Gamma_{k}(x)$, are presented in Fig. 7 using the data from Fig. 6.

The following result ensures that

[TABLE]

Theorem 3**.**

For any value of noise $\xi(x)\in[-\delta,\delta]$ the function $\Gamma_{k}(x)$ from (20) is a majorant for the objective function $g(x)$ satisfying (7) over the search region $D$, i.e., condition (22) holds.

Proof.

The theorem follows immediately from Theorem 2. ∎

Corollary 2**.**

Function $\Gamma_{k}(x)$ from (20) satisfies Lipschitz condition (1).

Proof.

It follows from considerations analogous to proof of Lemma 1 that functions $\gamma(x,x_{i})$, $1\leq i\leq k,$ satisfy Lipschitz condition (1). Since $\Gamma_{k}(x)$ is minimum of $k$ Lipschitz functions satisfying (1), this condition holds for $\Gamma_{k}(x)$, as well. ∎

In order to estimate the global maximizer of $g(x),x\in\widehat{\Omega},$ we propose a procedure that generalizes to the case of maximization of noisy functions the algorithm of Piyavskij (see [34]) proposed for Lipschitz optimization and working without any noise (this algorithm is called PA hereinafter). The PA works constructing a majorant $\Lambda_{k}(x)$ similar to $\Gamma_{k}(x)$ at each iteration using previously executed evaluations of the Lipschitz objective function $f(x)$. The absence of the noise allows the PA to improve $\Lambda_{k}(x)$ after each evaluation of $f(x)$ that is executed at a point

[TABLE]

In other words, due to the Lipschitz condition (1) it follows that

[TABLE]

and, as a consequence, after each new evaluation of $f(x)$ the majorant $\Lambda_{k}(x)$ becomes closer and closer to the objective function $f(x)$ allowing the PA to obtain a good approximation of the global maximum.

In contrast, in the noisy framework we deal with, a condition for $\Gamma_{k}(x)$ similar to (23) cannot be written down since the noise (see Fig. 6) can provoke a situation where

[TABLE]

In cases where the condition (24) holds, the evaluation of $g(x^{k+1})$ does not improve $\Gamma_{k}(x)$ and leads to the equivalence $\Gamma_{k+1}(x)=\Gamma_{k}(x)$. Let us recall that we have already faced similar situations before, when we got results of evaluation of $g(x)$ that did not allow us to enlarge the current safe region. Therefore, the remedy will be the same, i.e., execution of $\nu(x^{k+1})\leq\nu$ repeated evaluations of $g(x)$ at the point $x^{k+1}$ from (24) until either the number of possible repetitions $\nu$ will be reached or a value less than $\max_{x\in\widehat{\Omega}}\Gamma_{k}(x)$ will be produced. In the former case the search stops and in the latter one (similarly to (19)) the value

[TABLE]

is accepted and $\Gamma_{k+1}(x)\leq\Gamma_{k}(x)$ is built.

In practice, the operations of building $\Gamma_{k}(x)$ from the trial data $(x_{i},\check{g}(x_{i})),1\leq i\leq k,$ provided by Algorithm 1 are executed as follows. To each trial point $x_{i}$ we associate a value

[TABLE]

Then, it follows from the Lipschitz condition (1) that

[TABLE]

[TABLE]

and it is reached at the point

[TABLE]

where hereinafter the value $R_{i}$ is called characteristic of the interval $[x_{i-1},x_{i}]$. As a result, it follows that

[TABLE]

and, therefore, the point $x^{k+1}$ from (24) can be calculated as

[TABLE]

If the interval $[x_{t-1},x_{t}]$ is larger than the preset accuracy $\varepsilon$ (i.e., the stopping rule is not satisfied), the process of global maximization tries to update the majorant $\Gamma_{k}(x)$ by obtaining a value $\check{g}(x^{k+1})$ satisfying (25). In case this is not possible after $\nu$ evaluations of $g(x^{k+1})$, the algorithm stops. Otherwise, the interval $[x_{t-1},x_{t}]$ is subdivided in two subintervals $[x_{t-1},x^{k+1}]$ and $[x^{k+1},x_{t}]$, the number of points where $g(x)$ has been evaluated becomes $k+1$, the intervals having numbers larger than $t$ are renumbered and the two new intervals become $[x_{t-1},x_{t}]$ and $[x_{t},x_{t+1}]$, respectively. Then, the value $\check{g}(x^{k+1})$ is assigned to $z_{t}$ and all the values $z_{i}$ are renewed as follows

[TABLE]

This concludes construction of the new majorant $\Gamma_{k+1}(x)$ and the optimization process can be repeated. A detailed description of Algorithm 2 executing the global maximization over $\widehat{\Omega}$ is given in Appendix (the maximization is performed separately at each subregion of $\widehat{\Omega}$, but it can also be executed simultaneously over all disjoint subregions of the safe region, thus accelerating the search).

After satisfaction of one of the two stopping rules of Algorithm 2, a lot of information regarding the problem (1), (2), (3)–(5) is returned. First of all, the found safe region $\widehat{\Omega}$ is provided. Then, the largest found value $g^{*}_{k}$ and the corresponding point $x^{*}_{k}$, namely,

[TABLE]

can be taken as an estimate of the global maximum for the safe optimization problem (1), (17), (3)–(5). Then, the minorant $\Phi_{k}(x)$ and the majorant $\Gamma_{k}(x)$ for functions $f(x)$ and $g(x)$ are supplied, i.e.,

[TABLE]

where both $\Phi_{k}(x)$ and $\Gamma_{k}(x)$ are Lipschitz piece-wise linear functions (see Fig. 7). Their simple structure allows one to delimit easily areas of $\widehat{\Omega}$ where the global maximizers of $g(x)$ and $f(x)$ cannot be located. In fact, the set

[TABLE]

cannot contain the global maximizer of $g(x)$ and

[TABLE]

cannot contain the global maximizer of $f(x)$ (see Fig. 8 for illustration). The following lemma establishes that each new trial point produced by Algorithm 2 is chosen with the tentative to reduce the area $\widehat{\Omega}\setminus N^{k}_{g}$ containing the global maximizer of $g(x)$. This reduction becomes effective in case the condition (24) does not hold at the chosen point $x^{k+1}$.

Lemma 2**.**

All trial points produced by Algorithm 2 at any iteration number $k>1$ will belong to the area $\widehat{\Omega}\setminus N^{k}_{g}$.

Proof.

Due to its definition, the area $N^{k}_{g}$ contains points $x$ for which $\Gamma_{k}(x)<g^{*}_{k}$. It follows from (29) and (30) that the new trial point $x^{k+1}$ will be chosen by Algorithm 2 at the point corresponding to $\max_{x\in\widehat{\Omega}}\Gamma_{k}(x)$. The fact that this value is strictly larger than $g^{*}_{k}$ concludes the proof. ∎

The introduced Algorithm 2 works with the function $g(x)$ and looks for its global maximizer $x^{*}$ from (17). However, sometimes the interest of the person executing the optimization is related to the original function $f(x)$. Clearly, due to the noise, the global maxima of $g(x)$ and $f(x)$ can be different both in value and in location. The following theorem establishes conditions where Algorithm 2 can localize not only the area $\widehat{\Omega}\setminus N^{k}_{g}$ containing the global maximum of the function $g(x)$ but also can indicate a neighborhood of the global maximizer $X^{*}$ of $f(x)$ over $\widehat{\Omega}$.

Theorem 4**.**

Suppose that:

i. The function $f(x)$ satisfies the Lipschitz condition (1).

ii. The function $f(x)$ has the only global maximizer $X^{*}$ over $\widehat{\Omega}$ and for any local maximizer $x^{\prime}$ of $f(x)$ such that $x^{\prime}\neq X^{*},\,x^{\prime}\in\widehat{\Omega},$ it follows that

[TABLE]

where $\Delta>0$ is a fixed finite number and $\delta$, as usual, is the maximal level of the noise from (3).

iii. Parameters of Algorithm 2 are chosen as follows: $\varepsilon=0$, $\nu=\infty$.

Then during the work of Algorithm 2 there exists an iteration number $d^{*}$ such that for $k>d^{*}$ the only local maximizer of $f(x)$ located at the region $\widehat{\Omega}\setminus N^{k}_{g}$ is $X^{*}$.

Proof.

Supposition (iii) regarding parameters of Algorithm 2 means that the majorant $\Gamma_{k}(x)$ will be improved at each iteration at each new trial point $x^{k+1}$ because the number of allowed repeated evaluations of $g(x^{k+1})$ is $\nu=\infty$ and, therefore, the algorithm will not stop due to this stopping rule (see lines 48–49 of Algorithm 2). Then, since we have $\varepsilon=0$ in the second stopping rule, the algorithm will not stop due to this rule (see lines 43 and 70 of Algorithm 2), as well. Thus, supposition (iii) ensures that the sequence of trial points generated by Algorithm 2 is infinite and the majorant $\Gamma_{k}(x)$ is improved at each iteration.

Since $f(x)$ satisfies the Lipschitz condition (1) and due to the finiteness of $\Delta$ it follows from supposition (ii) and condition (33) that there exists a neighborhood $\omega(X^{*})$ of the point $X^{*}$ such that

[TABLE]

[TABLE]

and there are no other local maximizers of $f(x)$ belonging to $\omega(X^{*})$.

By taking into account the definition of $g(x)$ in (2) and (3), condition (34) can be re-written in the form

[TABLE]

meaning that for any point $x\in\omega(X^{*})$ and any local maximizer $x^{\prime}\neq X^{*},\,x^{\prime}\in\widehat{\Omega},$ it follows that the value of $f(x)$ decreased by the maximal level of the noise is larger than the value of $f(x^{\prime})$ increased by the highest level of the noise and, as a result, we have

[TABLE]

On the other hand, Algorithm 2 works in such a way that the majorant $\Gamma_{k}(x)$ is improved at each iteration and the point of its improvement corresponds to its maximum (see (29) and (30)). This means that $\max_{x\in\widehat{\Omega}}\Gamma_{k}(x)$ decreases with the increase of the iteration number $k$ (notice that this decrease is not strict since there can be situations where $\max_{x\in\widehat{\Omega}}\Gamma_{k}(x)$ is reached at several points). The decrease of $\max_{x\in\widehat{\Omega}}\Gamma_{k}(x)$ and the presence of a finite neighborhood $\omega(X^{*})$ where all points $x\in\omega(X^{*})$ satisfy (35) mean that there will exist an iteration number $\tilde{d}$ such that the new point $x^{\tilde{d}}$ will belong to $\omega(X^{*})$ and, therefore, the value of the function $g(x^{\tilde{d}})$ will satisfy

[TABLE]

for all local maximizers $x^{\prime}\neq X^{*},\,x^{\prime}\in\widehat{\Omega}$. Due to Lemma 2, new trials with the numbers $k>\tilde{d}$ will be executed at regions where $\Gamma_{k}(x)>g(x^{\tilde{d}})$. Therefore, since (36) holds and due to the already mentioned decrease of $\max_{x\in\widehat{\Omega}}\Gamma_{k}(x)$, after a finite number $\hat{d}$ of iterations of Algorithm 2 new trials will not be placed outside $\omega(X^{*})$. In fact, suppose that there exists a local maximizer $x^{\prime}\in\widehat{\Omega}$ that does not belong to $\omega(X^{*})$, at the $k$-th iteration $x^{\prime}\in[x_{i(k)-1},x_{i(k)}]$, and new trials are placed infinitely many times within this interval. It follows from (28) and footnote 1 (see section 2) that in this case

[TABLE]

and, therefore,

[TABLE]

As a result, for the characteristic $R_{i(k)}$ from (27) of this interval we have

[TABLE]

where the sign $\leq$ here is due to the rule (31). However, since $g(x^{\tilde{d}})>g(x^{\prime})$, the new trial cannot be placed within the interval $[x_{i(k)-1},x_{i(k)}]$ and, therefore, our supposition that new trial points will be placed in this interval infinitely many times is false.

Thus, by taking $d^{*}=\tilde{d}+\hat{d}$ and reminding that $\omega(X^{*})$ does not contain other local maximizers we conclude the proof. ∎

5 Numerical experiments

The experiments have been organized on three test cases, where the parameters of the $\delta$-Lipschitz framework have been set as follows: both $\delta$ and $h$ have been set at $10\%$ of the min-max range of $f(x)$, $\varepsilon=0.001$, $\nu=15$, and $\sigma=0.10$. For all test cases, initial safe points have been chosen randomly.

• •

Test case 1. This test case consists of four test problems. The safe expansion phase can lead to the identification of a unique safe region. Table 2 presents the structure of the test problems used in this test case. Fig. 9 shows for each test problem the final state of each phase in optimization processes (safe expansion – phase 1 and global maximization – phase 2).

• •

Test case 2. This test case consists of six test problems. The safe expansion phase can lead to the identification of the maximum of two disconnected safe regions. Table 2 presents the structure of the test problems used in this test case. Figs. 10–11 show for each test problem the final state of each phase in optimization processes (safe expansion – phase 1 and global maximization – phase 2).

• •

Test case 3. This test case consists of eight test problems. The safe expansion phase can lead to the identification of multiple disconnected safe regions. Table 3 presents the structure of the test problems used in this test case. Figs. 12–13 show for each test problem the final state of each phase in optimization processes (safe expansion – phase 1 and global maximization – phase 2).

Table 4 presents results of the numerical experiments and Figs. 9–13 illustrate them. Figures related to Phase 1 show initial safe points by green dots and results of evaluations by symbol ‘x’; they present also the constructed minorant $\Phi_{k}(x)$ and the found safe region. Figures related to Phase 2 show results of evaluations executed during the global maximization by blue symbols ‘x’. The best found point is shown by symbol ‘+’ and the global maximizer by symbol ‘*’. These figures show also the constructed majorant $\Gamma_{k}(x)$ together with the functions $\gamma(x,x_{i})$ from (21) used to build $\Gamma_{k}(x)$. The black dots above the majorant represent values that either have not improved the majorant or better values (see (31)) than these ones have been obtained during the maximization process.

6 Conclusions

The problem of safe global maximization of an expensive black-box function $f(x)$ satisfying the Lipschitz condition has been considered in this paper. The word expensive means here that each evaluation of the objective function $f(x)$ is a time consuming operation. With respect to the traditional formulations used in global optimization the problem under consideration here has two important distinctions:

1. i.

The first difficulty consists in the presence of the noise. As a result, instead of the function $f(x)$ a function $g(x)$ distorted by the noise is optimized. 2. ii.

The second difficulty is related to the notion of safe optimization. The term “safe” means that the objective function $g(x)$ during optimization should not violate a “safety” threshold, that in our case is a certain a priori given value $h$ in the maximization problem. Any new function evaluation must be performed at “safe points” only, namely, at points $y$ for which it is known that the objective function $g(y)\geq h$. Clearly, the main difficulty here consists in the fact that the used optimization algorithm should ensure that the safety constraint will be satisfied at a point $y$ before evaluation of $f(y)$ will be executed.

In our approach, the problem under investigation has been split in two parts (phases):

During the first phase, it is required to find an approximation $\widehat{\Omega}$ of the safe region $\Omega$ within the search domain $D$ by learning from the information received from evaluations of $g(x)$. 2. 2.

Then, during the second phase, it is necessary to find an approximation of the global maximum of $g(x)$ over $\widehat{\Omega}$.

A theoretical study of the problem has been performed and it has been shown that even though the objective function $f(x)$ satisfies the Lipschitz condition, traditional Lipschitz minorants and majorants cannot be used in this context due to the presence of the noise. In fact, a counterexample showing that the usage of simple Lipschitz ideas can lead to erroneous evaluations of $g(x)$ at unsafe points has been presented.

Then, a $\delta$-Lipschitz framework has been introduced and an algorithm determining the safe area within the search domain $D$ has been proposed. It has been theoretically proven that this method is able to construct a minorant ensuring that all executed evaluations of $g(x)$ performed by the method will be safe. It has been shown that the introduced procedure allows one to expand reliably the safe region from initial safe points and to obtain the desired approximation $\widehat{\Omega}$ of the safe region $\Omega$.

After that the second algorithm executing the safe global maximization over the found safe region $\widehat{\Omega}$ has been proposed. It has been theoretically investigated and conditions allowing one to make conclusions with respect to not only the noisy function $g(x)$ but also regarding the original function $f(x)$ have been established. It should be stressed that the introduced algorithm is able not only to provide lower and upper bounds for the global maxima of $g(x)$ and $f(x)$ but to restrict significantly the area of a possible location of the respective global maximizers, as well.

Finally, numerical experiments executed on a series of problems showing the reliability of the proposed procedures have been performed.

Compliance with Ethical Standards

Conflict of interest: The authors declare that they have no conflict of interest.

Ethical approval: This article does not contain any studies with human participants or animals performed by any of the authors.