Global Linear Convergence of Evolution Strategies on More than Smooth Strongly Convex Functions

Throughout the paper, we use the following notations. The set of natural numbers $\{1,2,\ldots,\}$ is denoted $\mathbb{N}$. Open, closed, and left open intervals on $\mathbb{R}$ are denoted by $(,)$, $[,]$, and $(,]$, respectively. The set of strictly positive real numbers is denoted by $\mathbb{R}_{>}$. The Euclidean norm on $\mathbb{R}^{d}$ is denoted by $\lVert~{}\rVert$. Open and closed balls with center $c$ and radius $r$ are denoted as $\mathcal{B}(c,r)=\{x\in\mathbb{R}^{d}:\lVert x-c\rVert<r\}$ and $\bar{\mathcal{B}}(c,r)=\{x\in\mathbb{R}^{d}:\lVert x-c\rVert\leqslant r\}$, respectively. Lebesgue measures on $\mathbb{R}$ and $\mathbb{R}^{d}$ are both denoted by the same symbol $\mu$. A multivariate normal distribution with mean $m$ and covariance matrix $\Sigma$ is denoted by $\mathcal{N}(m,\Sigma)$. Its probability measure and its induced probability density under Lebesgue measure are denoted by $\Phi(\cdot;m,\Sigma)$ and $\varphi(\cdot;m,\Sigma)$. The indicator function of a set or condition $C$ is denoted by $1\left\{\scriptstyle{C}\right\}$. We use Bachmann-Landau notations $f\in o(g)$, $O(g)$, $\omega(g)$, $\Omega(g)$ and $\Theta(g)$ to mean $\limsup_{\epsilon\to 0}f(\epsilon)/g(\epsilon)=0$, $\limsup_{\epsilon\to 0}f(\epsilon)/g(\epsilon)\leqslant C$ for some $C>0$, $\liminf_{\epsilon\to 0}f(\epsilon)/g(\epsilon)=\infty$, $\liminf_{\epsilon\to 0}f(\epsilon)/g(\epsilon)\geqslant C$ for some $C>0$, and $f\in O(g)$ and $f\in\Omega(g)$, respectively.

We analyze a generalized version of the (1+1)-ES with one-fifth success rule presented in Algorithm 1, which implements one of the oldest approaches to adapt the step-size in randomized optimization methods [51, 17, 54]. The specific implementation was proposed in [38]. At each iteration, a candidate solution $x_{t}$ is sampled. It is centered in the current incumbent $m_{t}$ and follows a multivariate normal distribution with mean vector $m_{t}$ and covariance matrix equal to $\sigma_{t}^{2}I_{d}$ where $I_{d}$ denotes the identity matrix. The candidate solution is accepted, that is $m_{t}$ becomes $x_{t}$, if and only if $x_{t}$ is better than $m_{t}$ (i.e. $f(x_{t})\leqslant f(m_{t})$). In this case, we say that the candidate solution is successful. The step-size $\sigma_{t}$ is adapted so as to maintain a probability of success to be approximately the target success probability denoted by $p_{\mathrm{target}}:=\frac{\log(1/\alpha_{\downarrow})}{\log(\alpha_{\uparrow}/\alpha_{\downarrow})}$. To do so, the step-size is increased by the increase factor $\alpha_{\uparrow}>1$ in case of success (which is an indication that the step-size is likely to be too small) and decreased by the decrease factor $\alpha_{\downarrow}<1$ otherwise. The covariance matrix $\Sigma_{t}$ of the sampling distribution of candidate solutions is adapted in the set $\mathcal{S}_{\kappa}$ of positive-definite symmetric matrices with determinant $\det(\Sigma)=1$ and condition number $\operatorname{Cond}(\Sigma)\leqslant\kappa$. We do not assume any specific update mechanism for $\Sigma$, but we assume that the update of $\Sigma$ is invariant to any strictly increasing transformation of $f$. We call such an update comparison-based (see Lines 7 and 11 of Algorithm 1). Then, our algorithm behaves exact-equally on $f$ and on $g\circ f$ for all strictly increasing functions $g:\mathbb{R}\to\mathbb{R}$ (i.e., $g(s)\lesseqgtr g(t)\Leftrightarrow s\lesseqgtr t$). This defines a class of comparison-based randomized algorithms and we denote it as (1+1)$\text{-ES}_{\kappa}$. For $\kappa=1$, it is simply denoted as (1+1)-ES.

Note that $\alpha_{\uparrow}$ and $\alpha_{\downarrow}$ are not meant to be tuned depending on the function properties. How to choose such constants for $\Sigma_{t}=I_{d}$ is well-known and is related to the so-called evolution window [52]. In practice, $\alpha_{\downarrow}=\alpha_{\uparrow}^{-1/4}$ is the most commonly used setting, which leads to $p_{\mathrm{target}}=1/5$. It has been shown to be close to optimal, which gives nearly optimal (linear) convergence rate on the sphere function [51, 17]. Hereunder we write $\theta=(m,\sigma,\Sigma)$ as the state of the algorithm, $\theta_{t}=(m_{t},\sigma_{t},\Sigma_{t})$ and the state-space is denoted by $\Theta$.

Figure 1 shows typical runs of the (1+1)-ES and a version of (1+1)$\text{-ES}_{\kappa}$ proposed in [6], which is known as the (1+1)-CMA-ES, on a $10$-dimensional ellipsoidal function with different condition numbers $\kappa_{f}$ of the Hessian. It is empirically observed that $\Sigma_{t}$ in the (1+1)-CMA-ES approaches the inverse Hessian $\nabla^{2}f(m_{t})$ of the objective function up to the scalar factor if the objective function is convex quadratic. The runtime of (1+1)-ES scales linearly with $\kappa_{f}$ (notice the logarithmic scale of the horizontal axis), while the runtime of the (1+1)-CMA-ES suffers only an additive penalty, roughly proportional to the logarithm of $\kappa_{f}$. Once the Hessian is well approximated by $\Sigma$ (up to a scalar factor), it approaches the global optimum geometrically at the same rate for different values of $\kappa_{f}$.

In our analysis, we do not assume any specific $\Sigma$ update mechanism, hence it does not necessarily behave as shown in Figure 1. Our analysis is therefore the worst case analysis (for the upper bound of the runtime) and the best case analysis (for the lower bound of the runtime) among the algorithms in (1+1)$\text{-ES}_{\kappa}$.

Given positive real numbers $a$ and $b$ satisfying $0\leqslant a<b\leqslant+\infty$, and a measurable objective function, let ${\mathcal{X}_{a}^{b}}$ be the set defined in Definition 2.6.

We pose two core assumptions on the objective functions under which we will derive an upper bound on the expected first hitting time of $[0,\epsilon]$ by $f_{\mu}(m_{t})$ (Theorem 4.6) provided $a\leqslant\epsilon\leqslant f_{\mu}(m_{0})\leqslant b$. First, we require to be able to embed and include balls of radius scaling with $f_{\mu}(m)$ into the sublevel sets of $f$. We do not require this to hold on the whole search space but, for a set $\mathcal{X}_{a}^{b}$.

1. A1

   We assume that $f$ is a measurable function and that there exists $\mathcal{X}_{a}^{b}$ such that for any $m\in\mathcal{X}_{a}^{b}$, there exist an open ball $\mathcal{B}_{\ell}$ with radius $C_{\ell}f_{\mu}(m)$ and a closed ball $\bar{\mathcal{B}}_{u}$ with radius $C_{u}f_{\mu}(m)$ such that it holds $\mathcal{B}_{\ell}\subseteq\{x\in\mathbb{R}^{d}\mid f_{\mu}(x)<f_{\mu}(m)\}$ and $\{x\in\mathbb{R}^{d}\mid f_{\mu}(x)\leqslant f_{\mu}(m)\}\subseteq\bar{\mathcal{B}}_{u}$.

We do not specify the center of those balls that may or may not be centered on an optimum of the function. We will see in Proposition 4.1 that this assumption allows to bound $p^{\mathrm{lower}}_{(a,b]}(\bar{\sigma})$ and $p^{\mathrm{upper}}_{(a,b]}(\bar{\sigma})$ by tractable functions of $\bar{\sigma}$ which will be essential for the analysis. The property is illustrated in Figure 2.

The second assumption requires that the functions are $p$-improvable for $p$ which is lower-bounded uniformly over ${\mathcal{X}_{a}^{b}}$.

1. A2

   Let $f$ be a measurable function, we assume that there exists ${\mathcal{X}_{a}^{b}}$ and there exists $p^{\mathrm{limit}}>p^{\mathrm{target}}$ such that for any $m\in\mathcal{X}_{a}^{b}$ and any $\Sigma\in\mathcal{S}_{\kappa}$, the objective function $f$ is $p$-improvable for some $p\geqslant p^{\mathrm{limit}}$, i.e.,

   (5)

   $$\liminf_{\bar{\sigma}\downarrow 0}p^{\mathrm{lower}}_{(a,b]}(\bar{\sigma})\geqslant p^{\mathrm{limit}}\enspace.$$

The property is illustrated in Figure 2. This assumption implies in particular for a continuous function that ${\mathcal{X}_{a}^{b}}$ does not contain any local optimum. This latter assumption is required to obtain global convergence [24, Theorem 2] even without any covariance matrix adaptation (i.e. with $\kappa=1$) and it can be intuitively understood: If we have a point which is $p$-improvable with $p<p_{\mathrm{target}}$ and which is not a local minimum of the function, then, starting with a small step-size, the success-based step-size control may keep decreasing the step-size at such a point and the (1+1)$\text{-ES}_{\kappa}$ will prematurely converge to a point that is not a local optimum.

If A1 is satisfied with balls centered at the optimum $x^{*}$ of the function $f$, then it is easy to see that for all $x\in{\mathcal{X}_{a}^{b}}$

If the balls are not centered at the optimum, we have the one-side inequality $\lVert x-x^{*}\rVert\leqslant 2C_{u}f_{\mu}(x)$. Hence, the expected first hitting time of $f_{\mu}(m_{t})$ to $[0,\epsilon]$ translates to an upper bound for the expected first hitting time of $\lVert m_{t}-x^{*}\rVert$ to $[0,2C_{u}\epsilon]$.

We remark that A1 and A2 satisfied for $a=0$ allow to include some non-differentiable functions with non-convex sublevel sets as illustrated in Figure 2.

We now give two examples of functions that satisfy A1 and A2, including function classes where linear convergence of numerical optimization algorithms are typically analyzed. The first class consists of quadratically bounded functions. It includes all strongly-convex functions with Lipschitz continuous gradient. It also includes some non-convex functions. The second class consists of positively homogeneous functions. The levelsets of a positively homogeneous function are all geometrically similar around $x^{*}$.

1. A3

   We assume that $f=g\circ h$ where $g$ is a strictly increasing function and $h$ is measurable, continuously differentiable with the unique critical point $x^{*}$, and quadratically bounded around $x^{*}$, i.e., for some $L_{u}\geqslant L_{\ell}>0$,

   (7)

   $\displaystyle(L_{\ell}/2)\lVert x-x^{*}\rVert^{2}\leqslant h(x)-h(x^{*})\leqslant(L_{u}/2)\lVert x-x^{*}\rVert^{2}\enspace.$

2. A4

   We assume that $f=g\circ h$ where $h$ is continuously differentiable and positively homogeneous with a unique optimum $x^{*}$, i.e., for some $\gamma>0$

   (8)

   $$h(x^{*}+\gamma x)=h(x^{*})+\gamma\left(h(x^{*}+x)-h(x^{*})\right)\enspace.$$

The following lemmas show that these assumptions imply A1 and A2. The proofs of the lemmas are presented in Section B.1 and Section B.2, respectively.

We start with the generic definition of the first hitting time of a stochastic process $\{X_{t}:t\geqslant 0\}$, defined as follows.

The first hitting time is the number of iterations that the stochastic process requires to reach the target level $\beta<\beta_{0}$ for the first time. In our situation, $X_{t}=\lVert m_{t}-x^{*}\rVert$ measures the distance from the current solution $m_{t}$ to the target point $x^{*}$ (typically, global or local optimal point) after $t$ iterations. Then, $\beta=\epsilon>0$ defines the target accuracy and $T_{\epsilon}^{X}$ is the runtime of the algorithm until it finds an $\epsilon$-neighborhood $\mathcal{B}(x^{*},\epsilon)$. The first hitting time $T_{\epsilon}^{X}$ is a random variable as $m_{t}$ is a random variable. In this paper, we focus on the expected first hitting time $\mathbb{E}[T_{\epsilon}^{X}]$. We want to derive lower and upper bounds on this expected hitting time that relate to the linear convergence of $X_{t}$ towards $x^{*}$. Such bounds take the following form: There exist $C_{T},\tilde{C}_{T}\in\mathbb{R}$ and $C_{R}>0$, $\tilde{C}_{R}>0$ such that for any $0<\epsilon\leqslant\beta_{0}$

That is, the time to reach the target accuracy scales logarithmically with the ratio between the initial accuracy $\lVert m_{0}-x^{*}\rVert$ and the target accuracy $\epsilon$. The first pair of constants, $C_{T}$ and $\tilde{C}_{T}$, capture the transient time, which is the time that adaptive algorithms typically spend for adaptation. The second pair of constants, $C_{R}$ and $\tilde{C}_{R}$, reflect the speed of convergence (logarithmic convergence rate). Intuitively, assuming that $C_{R}$ and $\tilde{C}_{R}$ are close, the distance to the optimum decreases in each step at a rate of approximately $\exp(-C_{R})\approx\exp(-\tilde{C}_{R})$. While upper-bounds inform us about the (linear) convergence, the lower-bound helps understanding whether the upper bounds are tight.

Alternatively, linear convergence can be defined as the following property: there exits $C>0$ such that

When we have an equality in the previous statement, we say that $\exp(-C)$ is the convergence rate.

Figure 1 (right plot) visualizes three different runs of the (1+1)-ES on a function with spherical level sets with different initial step-sizes. First of all, we clearly observe linear convergence. The first hitting time of $\mathcal{B}(x^{*},\epsilon)$ scales linearly with $\log(1/\epsilon)$ for a sufficiently small $\epsilon>0$. Second, its convergence speed is independent of the initial condition. Therefore, we expect to have universal constants $C_{R}$ and $\tilde{C}_{R}$ independent of the initial state. Last, depending on the initial step-size, the transient time can vary. If the initial step-size is too large or too small, it does not produce progress in terms of $\lVert m_{t}-x^{*}\rVert$ until the step-size is well adapted. Therefore, $C_{T}$ and $\tilde{C}_{T}$ depend on the initial condition, with a logarithmic dependency on the initial multiplicative mismatch.

We are going to use drift analysis that consists in deducing properties of a sequence $\{X_{t}:t\geqslant 0\}$ (adapted to a natural filtration $\{\mathcal{F}_{t}:t\geqslant 0\}$) from its drift defined as $\mathbb{E}[X_{t+1}\mid\mathcal{F}_{t}]-X_{t}$ [29]. Drift analysis has been widely used to analyze hitting times of evolutionary algorithms defined on discrete search spaces (mainly on binary search spaces) [32, 33, 11, 46, 20, 19]. Though they were developed mainly for finite search spaces, the drift theorems can naturally be generalized to continuous domains [42, 44]. Indeed, Jägersküpper’s work [34, 36, 37] is based on the same idea, while the link to the drift analysis was implicit.

Since many drift conditions have been developed for analyzing algorithms on discrete domains, the domain of $X_{t}$ is often implicitly assumed to be bounded. However, this assumption is violated in our situation, where we will use $X_{t}=\log\big{(}f_{\mu}(m_{t})\big{)}$ as the quality measure, which takes values in $\mathbb{R}\cup\{-\infty\}$, and is meant to approach $-\infty$. We refer to [4] for additional details. In general, translating expected progress requires bounding the tail of the progress distribution, as formalized in [29].

To control the tails of the drift distribution, we construct a stochastic process $\{Y_{t}:t\geqslant 0\}$ iteratively as follows: $Y_{0}=X_{0}$ and

for some $A\geqslant B>0$ and $\beta<\beta_{0}$ with $X_{0}=\beta_{0}$. It clips $X_{t+1}-X_{t}$ to some constant $-A$ ($A>0$) from below. We introduce the indicator $1\left\{\scriptstyle{T_{\beta}^{X}>t}\right\}$ for a technical reason. The process disregards progress larger than $A$, and it fixes the progress of the step that hits the target set to $B$. It is formalized in the following theorem, which is our main mathematical tool to derive an upper bound of the expected first hitting time of (1+1)$\text{-ES}_{\kappa}$ in the form of (9).

This theorem can be intuitively understood: we assume for the sake of simplicity a process $X_{t}$ such that $X_{t+1}\geqslant X_{t}-A$. Then (12) states that the process progresses in expectation by at least $-B$. The theorem concludes that the expected time needed to reach a value smaller than $\beta$ when started in $\beta_{0}$ equals to $(\beta_{0}-\beta)/B$ (what we would get for a deterministic algorithm) plus $A/B$. This last term is due to the stochastic nature of the algorithm. It is minimized if $A$ is as close as possible to $B$, which corresponds to a highly concentrated process.

Jägersküpper [36, Theorem 2] established a general lower bound of the expected first hitting time of the (1+1)-ES. We borrow the same idea to prove the following general theorem for a lower bound of the expected first hitting time, which generalizes [37, Lemma 12]. See Theorem 2.3 in [4] for its proof.

In the sequel, we will analyze the process $\{\theta_{t}:t\geqslant 0\}$ where $\theta_{t}=(m_{t},\sigma_{t},\Sigma_{t})\in\mathbb{R}^{n}\times\mathbb{R}_{>}\times\mathcal{S}_{\kappa}$ generated by the (1+1)$\text{-ES}_{\kappa}$ algorithm. We assume from now on that the optimized objective function $f$ is measurable with respect to the Borel $\sigma$-algebra. We equip the state-space $\mathcal{X}=\mathbb{R}^{n}\times\mathbb{R}_{>}\times\mathcal{S}_{\kappa}$ with its Borel $\sigma$-algebra denoted $\mathcal{B}(\mathcal{X})$.

We present two preliminary results. In Assumption A1, we assume that for $m\in{\mathcal{X}_{a}^{b}}$, we can include a ball of radius $C_{\ell}f_{\mu}(m)$ into the sublevel set $S_{0}(m)$ and embed $S_{0}(m)$ into a ball of radius $C_{u}f_{\mu}(m)$. This allows us to upper bound and lower bound the probability of success for all $m\in{\mathcal{X}_{a}^{b}}$, for all $\Sigma\in\mathcal{S}_{\kappa}$, by the probability to sample inside of balls of radius $C_{u}f_{\mu}(m)$ and $C_{\ell}f_{\mu}(m)$ with appropriate center. From this we can upper-bound $p^{\mathrm{upper}}_{(a,b]}(\bar{\sigma})$ by a function of $\bar{\sigma}$. Similarly we can lower-bound $p^{\mathrm{lower}}_{(a,b]}(\bar{\sigma})$ by a function of $\bar{\sigma}$. The corresponding proof is given in Section B.3.