# Imaginary Zeroth-Order Optimization

Wouter Jongeneel111The author is with the Risk Analytics and Optimization
Chair, École Polytechnique Fédérale de Lausanne (EPFL) and is supported by the
Swiss National Science Foundation under the NCCR Automation, grant agreement
51NF40_180545. The author is grateful to RAO colleagues and ICCOPT 2022
participants, their comments greatly improved the work. Contact:
<EMAIL_ADDRESS>wjongeneel.nl

(upload: December 14, 2021, last update: August 14, 2022)

###### Abstract

Zeroth-order optimization methods are developed to overcome the practical
hurdle of having knowledge of explicit derivatives. Instead, these schemes
work with merely access to noisy functions evaluations. One of the predominant
approaches is to mimic first-order methods by means of some gradient
estimator. The theoretical limitations are well-understood, yet, as most of
these methods rely on finite-differencing for shrinking differences, numerical
cancellation can be catastrophic. The numerical community developed an
efficient method to overcome this by passing to the complex domain. This
approach has been recently adopted by the optimization community and in this
work we analyze the practically relevant setting of dealing with computational
noise. To exemplify the possibilities we focus on the strongly-convex
optimization setting and provide a variety of non-asymptotic results,
corroborated by numerical experiments, and end with local non-convex
optimization.

Keywords—zeroth-order optimization, derivative-free optimization, complex-step
derivative, gradient estimation, numerical optimization.
AMS Subject Classification (2020)—65D25, 65G50, 65K05, 65Y04, 65Y20, 90C56.

##  1 Introduction

“La voie la plus courte et la meilleure entre deux vérités du domaine réel
passe souvent par le domaine imaginaire.”—J. Hadamard222See
http://homepage.math.uiowa.edu/~jorgen/hadamardquotesource.html.

From the Fourier transformation, quantum mechanics to the Nyquist stability
criterion, the complex numbers grew out to be quintessential mathematical
machinery.

Building upon the work by [KW52, LM67, NY83, ST98, FKM04, NS17], it is shown
in [JYK21] that (randomized) zeroth-order optimization also benefits from
passing to the complex domain as one can derive an inherently numerically
stable method, which is in sharp contrast to common finite-difference methods.
This work departs from [JYK21] by introducing an indispensable layer of
realism; noise.

We are interested in numerically solving optimization problems of the form

$\operatorname*{\text{minimize}}_{x\in\mathcal{X}}\quad f(x),$

where $f:\mathcal{D}\to\mathbb{R}$ is a smooth objective function defined on
an open set $\mathcal{D}\subseteq\mathbb{R}^{n}$, and
$\mathcal{X}\subseteq\mathcal{D}$ is a non-empty closed feasible set.
Optimizers, which based on the context could be globally or locally optimal,
are denoted by $x^{\star}$. We extend [JYK21] and assume that the objective
function $f$ can only be accessed through a zeroth-order oracle that outputs
corrupted function evaluations at prescribed test points, that is, with noise.
As we only have access to such a zeroth-order oracle, our work belongs to the
field of zeroth-order optimization, derivative-free optimization or more
generally black-box optimization [CSV09, AH17a].

We start by highlighting two important assumptions made throughout this work.

######  Assumption 1.1 (Smoothness).

The objective function $f$ is real-analytic over
$\mathcal{D}\subseteq\mathbb{R}^{n}$, denoted $f\in C^{\omega}(\mathcal{D})$.

Recall, a function is real-analytic when it can be locally expressed by a
convergent power series, which is stronger than smoothness, i.e.,
$C^{\omega}(\mathcal{D})\subset C^{\infty}(\mathcal{D})$. A complex-analytic
function is called holomorphic333More formally, a complex differentiable
function is called holomorphic, but as it turns out, complex differentiability
coincides with complex analyticity [Kra00].. With few exceptions [AMA05],
Assumption 1.1 does not appear often explicitly in the optimization
literature. However, by means of the results in [Pol86] it does appear
indirectly in for example the context of reinforcement-learning [Faz+18,
Mal+19]. As in [JYK21], Assumption 1.1 is again mainly there to allow for the
next assumption. As will be explained below, having access to $\Im(f(z))$ for
some $z\in\mathbb{C}^{n}$ is at the core of the approach. In contrast to
[JYK21] we allow for the presence of (computational) noise.

######  Assumption 1.2 (Stochastic complex oracle).

Consider some unknown function $f\in C^{\omega}(\mathcal{D})$ which admits a
holomorphic extension to $\Omega\subseteq\mathbb{C}^{n}$. We assume to have
access to an oracle which can output $\Re(f(z))+\xi$ and $\Im(f(z))+\xi$ for
any $z\in\Omega$ with $\xi$ a zero-mean random variable supported on
$\Xi\subseteq\mathbb{R}$ with $\mathbb{E}[\xi^{2}]\leq\sigma_{\xi}$ for some
$\sigma_{\xi}>0$.

Assumption 1.2 is particularly important in the simulation-based context. As
there the evaluation of $f(z)$ might pertain to millions of floating-point
operations, chopping and round-off errors are easily introduced. The set
$\Omega$ will be specified later on. We will make no further assumptions
regarding the distribution of $\xi$.

###  1.1 Related work

Arguably the first algorithm that uses noisy finite-differences to approximate
gradient algorithms is the Kiefer-Wolfowitz algorithm [KW52], [KC78, Section
2.3.5]. [NY83] contributed the first single-point gradient estimator and
perhaps more importantly, the need for lower bounds. A large fraction of the
work on zeroth-order optimization entails mimicking first-order algorithms via
some approximation of the gradient. These types of algorithms are generally
scalable444See however the discussion in [Sch22] to put this in the correct
perspective., easy to implement and as they mimic first-order methods, they
usually come with guarantees. A common gradient estimator is of the form

$\widehat{g}_{\delta}(x)=\frac{n}{\delta}\left(f(x+\delta u)-f(x)\right)u,$
(1.1)

for some choice of the smoothing parameter $\delta>0$ and some appropriately
chosen random variable $u$. The $\delta$ is sometimes referred to as the
exploration parameter. See that (1.1) requires two function evaluations, as
such we speak of a multi-point method. Using estimators of the form (1.1) was
popularized in the bandit-context [FKM04], although for a single-point
estimator, and relates largely to work on stochastic approximation algorithms
[KY03, Spa05] and to some extent to inexact/biased first-order methods [LT93,
d’A08, DGN14, AS21], [TSAK21, Section 4].

Compared to first-order methods, zeroth-order methods are commonly $O(n)$
times slower in the deterministic setting [NS17]. When noise is involved the
balance between bias and variance requires a more careful selection of the
smoothing parameter $\delta$. Let $K\in\mathbb{N}$ denote the length of the
sequence $x_{1},x_{2},\dots,x_{K}$ designed with the aim of converging in some
sense to (some) $\operatorname*{argmin}_{x\in\mathcal{X}}f(x)$. Let
$\bar{x}_{K}=K^{-1}\sum^{K}_{k=1}x_{k}$ be a uniformly-averaged iterate, we
will be mostly interested in quantifying how fast the optimization error

$\mathbb{E}[f(\bar{x}_{K})-f(x^{\star})]$

decays. Here, the expectation $\mathbb{E}[\cdot]$ is over the oracle noise and
the deliberate randomization within the proposed algorithms. That is, we can
define an abstract probability space $(\Omega,\mathcal{F},\mathbb{P})$ and
define $\mathbb{E}[\cdot]$ to be the expectation with respect to $\mathbb{P}$.
In [JNR12] the authors consider $\tau$-strongly convex functions with
$L$-Lipschitz gradients and show that the expected optimization error decays
like $\Omega(\sqrt{n/K})$ when using noisy single-point oracles. [Sha13] shows
that in the quadratic case the result can be improved. If the objective $f$ is
$r$-times continuously differentiable, Chen shows that a rate of the order
$O(K^{-(r-1)/2r})$ is optimal [Che88]. If $f$ is strongly convex the optimal
rate becomes $O(K^{-(r-1)/r})$ [PT90]. See also [RSS12] for more on optimal
rates in the stochastic setting.

In [Duc+15] the authors show the information-theoretic optimality of multi-
point (two-point) methods, yet, in [JYK21] the authors show the numerical
superiority of single-point schemes. This work sets out to show to what extent
this observation prevails when noise is present. As highlighted throughout the
recent survey article by [LMW19], it is not clear if there is a single-point
method which is as fast as multi-point methods. This observation motivates
[Zha+22] to use some form of memory such that their estimator only demands a
single new point each call. Nevertheless, in the end their method is
reminiscent of a multi-point method. Another recent work observes how the
continuous-time notion of extremum seeking can be translated to a zeroth-order
optimization algorithm [CTL22]. Their method turns out to be a combination of
the aforementioned residual-feedback and momentum and achieves an optimization
error of the order $O(n/K^{2/3})$, for a restrictive class of problems and a
deterministic oracle. To the best of our knowledge, we will provide the first
real single-point method which is capable of achieving an optimal rate.

We focus on one particular approach to zeroth-order optimization. Different
and successful lines of attack relate to model-based (trust-region) [CSV09],
Bayesian [Moc12] and more broadly black-box optimization [AH17a].

#### Contribution

We show that catastrophic numerical cancellation errors are also inevitable in
the widely used noisy multi-point case. We will show that this non-
deterministic setting also benefits from the imaginary gradient estimator as
proposed in [JYK21]. Using this single-point estimator and building upon
[HRB08, APT20], we provide the non-asymptotic analysis for a variety of
algorithms. Specifically, we consider for strongly convex functions the
unconstrained, constrained, online and quadratic cases. In the last setting we
can show that the algorithm is rate-optimal. To comply with zeroth-order
knowledge we also propose an estimation scheme for the strong-convexity
parameter. As an outlook we provide a local result in the nonconvex case and
showcase PDE-constrained optimization as an area of application. Besides, we
generalize some results from [JYK21] and we hope that a secondary contribution
of this work is to bring numerical intricacies further to the attention.

#### Structure

We start in Section 2 by detailing numerical problems in zeroth-order
optimization. In Section 3 we highlight the imaginary gradient estimator as
proposed in [JYK21] to overcome the aforementioned obstacles. Section 4 and
Section 5 provide all algorithms, corresponding convergence rates and a few
numerical experiments. Section 6.2 briefly comments on merely smooth non-
analytic functions and we conclude the work in Section 6.4. Some auxiliary
results can be found in the appendix.

#### Notation

The real and imaginary parts of a complex number $z=a+ib\in\mathbb{C}$ are
denoted by $\Re(z)=a$ and $\Im(z)=b$, while
$\mathbb{B}^{n}=\\{x\in\mathbb{R}^{n}:\|x\|_{2}\leq 1\\}$ is the Euclidean
$n$-ball and $\mathbb{S}^{n-1}=\partial\mathbb{B}^{n}$ denotes the Euclidean
$(n-1)$-sphere. Let $\mathsf{Y}\subset\mathbb{R}^{n}$ be a Borel measurable
set such that $\partial\mathsf{Y}$ is an orientable compact differentiable
manifold. We write $y\sim\mathsf{Y}$ to declare that $y$ is a random vector
following the uniform distribution on $\mathsf{Y}$, and for any Borel
measurable function $g:\mathsf{Y}\subset\mathbb{R}^{n}\rightarrow\mathbb{R}$
we denote by

$\mathbb{E}_{y\sim\mathsf{Y}}[g(y)]=\frac{1}{\mathrm{vol}(\mathsf{Y})}\int_{\mathsf{Y}}g(y)\mathrm{d}V(y)$

the expected value of $g(y)$, where $\mathrm{d}V$ represents the Borel measure
induced by the volume form on $\mathsf{Y}$, and $\mathrm{vol}(\mathsf{Y})$
represents the volume of $\mathsf{Y}$. The set of all $r$ times continuously
differentiable real-valued functions on the open set
$\mathcal{D}\subseteq\mathbb{R}^{n}$ is denoted by $C^{r}(\mathcal{D})$. Non-
negative constants are denoted by $C_{1},C_{2},\dots,C_{y},C_{z}$. Their
values can change from line to line. Regarding complexity notation,
$\Omega(\cdot)$, $\Theta(\cdot)$ and $O(\cdot)$ have their usual meaning with
$\widetilde{O}(\cdot)$ hiding logarithmic factors. The proof contain explicit
errors, whenever possible.

Using the notation from [Nes03] a function $f$ is said to be
$C^{k,r}_{L}(\mathcal{D})$-smooth when $f$ is $k$ times continuously
differentiable with additionally having its $r^{\mathrm{th}}$-derivative being
$L$-Lipschitz over some open set $\mathcal{D}\subseteq\mathbb{R}^{n}$. Here,
$k$ is an element of $\mathbb{N}_{\geq 0}\cup\\{\infty\\}\cup\\{\omega\\}$.
That is, if $f\in C^{1,1}_{L_{1}(f)}(\mathcal{D})$, then, $f$ has a Lipschitz
gradient, i.e.,

$\|\nabla f(x)-\nabla f(y)\|_{2}\leq L_{1}(f)\|x-y\|_{2},\quad\forall
x,y\in\mathcal{D}.$ (1.2)

Similarly, if $f\in C^{2,2}_{L_{2}(f)}(\mathcal{D})$, then, $f$ has a
Lipschitz Hessian, i.e.,

$\|\nabla^{2}f(x)-\nabla^{2}f(y)\|_{2}\leq L_{2}(f)\|x-y\|_{2}\quad\forall
x,y\in\mathcal{D}.$ (1.3)

Instead of the $\ell_{2}$-norm $\|\cdot\|_{2}$ one can generalize the above to
any norm $\|\cdot\|$ and its dual $\|\cdot\|_{\star}$. Note that when $f\in
C^{\omega}(\mathcal{D})$, then, the existence of $L_{1}(f)$ and $L_{2}(f)$ for
$f$ restricted to compact subsets of $\mathcal{D}$ is trivial. Yet, to aid the
reader, we will always indicate when we work with these constants.

##  2 Numerical stability in zeroth-order optimization

Multi-point finite-difference estimators dominate the zeroth-optimization
literature, e.g., see [HL14, Duc+15, NS17, Gas+17, Sha17, APT20, LLZ21, NG21]
or the recent survey articles [LMW19, Liu+20]. The motivation largely follows
from the observation that the initial single-point schemes as proposed in
[NY83, FKM04] have an unbounded variance, even when the function evaluations
come without noise. The multi-point schemes avoid this by constructing
estimators akin to numerical directional derivatives [ADX10, NS17].

Nevertheless, as pointed out in [JYK21], multi-point schemes do suffer from
catastrophic numerical cancellation. See also [Shi+21] for an extensive
numerical study on the numerical performance of finite-difference methods in
the context of optimization.

###  2.1 Numerical cancellation

The smallest $\epsilon_{M}\in\mathbb{Q}_{>0}$ such that on a particular
machine $1+\epsilon_{M}>1$ is called the machine precision. Nowadays, the
number $\epsilon_{M}$ is commonly of the order $10^{-16}$, which is the number
we will use. So in general, for a continuous function
$f:\mathbb{R}\to\mathbb{R}$, when $x,y\in\mathbb{R}$ are chosen such that
$f(x)-f(y)\leq\epsilon_{M}$ the numerical evaluation of $f(x)-f(y)$ can be
problematic. Now for zeroth-order gradient estimators, given some $f\in
C^{r+1}(\mathbb{R})$ with $r\geq 1$, then in the approximation

$\partial_{x}f(x)=\frac{f(x+\delta)-f(x)}{\delta}+O(\delta)$ (2.1)

one cannot make $\delta>0$ arbitrarily small and expect to recover
$\partial_{x}f(x)$. For a sufficiently small $\delta$ the evaluations
$f(x+\delta)$ and $f(x)$ will be numerically indistinguishable and
cancellation errors appear, see [Ove01, Chapter 11]. Running into these
machine-precision problems is inherent to finite-difference (multi-point)
optimization methods as one looks for (at) the flattest part of $f$.

A celebrated work-around in the numerical community is the so-called complex-
step method. This approach was introduced in [LM67] with the first concrete
complex-step approach appearing in [ST98] and with later elaborations to
higher-order derivatives, matrices and Lie groups in [MSA03, AMH10, ASM15,
Abr+18, CWF20]. In short, via the Cauchy-Riemann equations one can show that
for a holomorphic function $f$, one has

$\partial_{x}f(x)=\frac{\Im(f(x+i\delta))}{\delta}+O(\delta^{2}).$ (2.2)

Not only is numerical cancellation impossible, the error term improved
compared to (2.1). This approach recently surfaced in the optimization
community [NS18, HS21] with the first complete deterministic non-asymptotic
analysis appearing in [JYK21]. The first applications of the complex-step
derivative to Reinforcement Learning appeared in [WS21, WZS21]. Of course, as
complex arithmetic is more expensive than real arithmetic, numerical stability
does not come for free555For example, to compute the multiplication of 2
complex numbers $(a_{1}+ib_{1})(a_{2}+ib_{2})$ one needs 3 real
multiplications ,i.e., $(a_{1}+b_{1})(a_{2}+b_{2})$, $a_{1}a_{2}$ and
$b_{1}b_{2}$, see also [AL81]..

To visualize the power of the complex-step approach we provide a short
example.

######  Example 2.1 (Numerical estimator stability).

We showcase the forward-difference ($\mathsf{fd}$), central-difference
($\mathsf{cd}$) and complex-step ($\mathsf{cs}$) for $f(x)=\log(x)$ at $x=1$,
that is, we compare

$\displaystyle f_{\mathsf{fd}}(x,\delta)$
$\displaystyle=\frac{f(x+\delta)-f(x)}{\delta},$ $\displaystyle
f_{\mathsf{cd}}(x,\delta)$
$\displaystyle=\frac{f(x+\delta)-f(x-\delta)}{2\delta},$ $\displaystyle
f_{\mathsf{cs}}(x,\delta)$
$\displaystyle=\frac{\Im\left(f(x+i\delta)\right)}{\delta}$

for $\delta\downarrow 0$, see Figure 1(a).

(a) Example 2.1, gradient estimator comparisons. See also [JYK21, Example
2.6].

(b) Suboptimality gap for Example 4.8, Theorem 4.7 (CS) vs. [APT20, Theorem
5.1] ($\beta$).

(c) Suboptimality gap for Example 4.11, Algorithm 4.1 (b) combined with the
estimation scheme (4.16).

Figure 2.1: Numerical experiments. Each figure displays all available data.

Only the complex-step estimator can reach machine precision, yet the other two
methods are used frequently in zeroth-order optimization under the assumption
that one can select $\delta$ arbitrarily close to $0$. As such, these methods
leave something to be desired, numerically.

At last we elaborate on Example 2.1 and visualize the imaginary lifting of
$f(x)$. That is, for $f(x)=x^{p}$, with $x\in\mathbb{R}$ and $p\in\mathbb{N}$,
we show $\Im(f(x+iy)/y)$. Indeed, for sufficiently small $y$ we see in Figure
2.2 that this number converges to $\partial_{x}f(x)$ for $p\downarrow 2$666See
http://wjongeneel.nl/ZO.gif for an animated version of Figure 2.2..

(a) $p=50$

(b) $p=25$

(c) $p=10$

(d) $p=2$

Figure 2.2: Further visualizations of (2.2), i.e., $\partial_{x}f(x)\approx\Im
f(x+iy)/y$.

##  3 Imaginary gradient estimation

In this section we summarize the main tool as set forth by [JYK21]. Motivated
by Example 2.1, we consider the imaginary $\delta$-smoothed version of $f$ as
proposed in [JYK21], that is

$f_{\delta}(x)=\mathbb{E}_{v\sim\mathbb{B}^{n}}\left[\Re\big{(}f(x+i\delta
v)\big{)}\right].$ (3.1)

Here, the parameter $\delta\in\mathbb{R}_{>0}$ is the tuneable smoothing
parameter and relates to the radius of the ball we average over. As mentioned
before, the offset $\delta v$ in (3.1) relates to exploration777This notion of
exploration could be a benefit of these randomized approaches [Sch22]., due to
our limited amount of information on the objective, each direction is
potentially worthwhile exploring and as such we consider a perfectly symmetric
shape; the ball $\mathbb{B}^{n}$. See [HL14] and Lemma 3.8 for comments and
results beyond $\mathbb{B}^{n}$.

To make sure $f_{\delta}$ is well-defined, $f(x+i\delta v)$ needs to be well-
defined and as such we assume the following.

######  Assumption 3.1 (Holomorphic extension).

The function $f:\mathcal{D}\subseteq\mathbb{R}^{n}\to\mathbb{R}$ is real-
analytic over the open set $\mathcal{D}$ and admits a holomorphic extension to
$\mathcal{D}\times
i\cdot(-\bar{\delta},\bar{\delta})^{n}\subset\mathbb{C}^{n}$ for some
$\bar{\delta}\in(0,1)$.

See [JYK21, Section 2.1] for more on the existence of such an extension. Note,
the interval $(0,1)$ is merely a convenient choice for the exposition.

Next we highlight the approximation quality of $f_{\delta}$.

######  Lemma 3.2 (Approximation quality of the complex-step function [JYK21,
Proposition 3.2]).

Let $f\in C^{\omega,1}_{L_{1}(f)}(\mathcal{D})$ with $L_{1}(f)\geq 0$ satisfy
Assumption 3.1 for some $\bar{\delta}\in(0,1)$. Then, for $f_{\delta}$ as in
(3.1) and any fixed $x\in\mathcal{D}$ and $\kappa\in(0,1)$ there exists some
constant $C_{0}\geq 0$, vanishing with $L_{1}(f)$, such that

$\displaystyle\left|f_{\delta}(x)-f(x)\right|$
$\displaystyle\leq{C_{0}\delta^{2}}\quad\forall\,\delta\in(0,\kappa\bar{\delta}].$
(3.2)

It is imperative to remark that convexity of $f$ does not always carry over to
$f_{\delta}$, e.g., see [JYK21, Example 3.6].

Now we state one of the key contributions of [JYK21], which is the integral
representation of $\nabla f_{\delta}$. This result is the complex-step version
of the approach as proposed in [NY83, Section 9.3] and popularized by [FKM04,
Lemma 1].

######  Lemma 3.3 (The gradient of the complex-step function [JYK21,
Proposition 3.3]).

Let $f\in C^{\omega}(\mathcal{D})$ satisfy Assumption 3.1 for some
$\bar{\delta}\in(0,1)$, then, $f_{\delta}$ as in (3.1) is differentiable and
for any $x\in\mathcal{D}$ we have for any $\delta\in(0,\bar{\delta})$

$\nabla
f_{\delta}(x)=\frac{n}{\delta}\cdot\mathbb{E}_{u\sim\mathbb{S}^{n-1}}\left[\Im\left(f(x+i\delta
u)\right)u\right].$ (3.3)

The following result allows for showing consistency, i.e.,
$\lim_{\delta\downarrow 0}\nabla f_{\delta}(x)=\nabla f(x)$.

######  Lemma 3.4 (Integration over the $(n-1)$-sphere).

Given any $x\in\mathbb{R}^{n}$, then

$\frac{n}{\mathrm{vol}(\mathbb{S}^{n-1})}\cdot\int_{\mathbb{S}^{n-1}}\langle
x,u\rangle u\mathrm{d}V(u)=x.$ (3.4)

Although this result is well-known, for completeness we also provide the
proof.

###### Proof.

First, rewrite (3.4) as $n\cdot\int_{\mathbb{S}^{n-1}}uu^{\mathsf{T}}du\,x$
and recall that $uu^{\mathsf{T}}x=\langle x,u\rangle u$. Now we would like to
show that
$n\cdot\int_{\mathbb{S}^{n-1}}uu^{\mathsf{T}}du=\mathrm{vol}(\mathbb{S}^{n-1})\cdot
I_{n}$. To that end, use the geometric tracing identity
$n\cdot\int_{\mathbb{S}^{n-1}}\langle Xu,u\rangle
du=\mathrm{Tr}(X)\cdot\mathrm{vol}(\mathbb{S}^{n-1})$ [GHL04, Lemma 3.100],
differentiating both sides with respect to $X$ yields
$n\cdot\int_{\mathbb{S}^{n-1}}uu^{\mathsf{T}}du=\mathrm{vol}(\mathbb{S}^{n-1})\cdot
I_{n}$ indeed, which concludes the proof. ∎

Since $f$ is real-analytic, the directional derivative at $x\in\mathcal{D}$ in
the direction $u\in\mathbb{S}^{n-1}$ is well-defined and given by
$\langle\nabla f(x),u\rangle$. Then, observe from (3.3) that the approximation
is asymptotically consistent, that is, by appealing to the dominated
convergence theorem we have

$\lim_{\delta\downarrow 0}\nabla
f_{\delta}(x)=\frac{n}{\mathrm{vol}(\mathbb{S}^{n-1})}\cdot\int_{\mathbb{S}^{n-1}}\langle\nabla
f(x),u\rangle u\mathrm{d}V(u)\overset{\eqref{equ:int:sphere}}{=}\nabla f(x).$
(3.5)

Showing consistency of this type, albeit for the estimator, was one of the key
observations in [ADX10, NS17] to reduce gradient estimator variance. Such an
observation does not hold for other known single-point estimators cf. [FKM04,
Section 1.1].

Lemma 3.3 provides us immediately with a (noisy) single-point estimator of
$\nabla f_{\delta}(x)$, namely

$g_{\delta}(x)=\frac{n}{\delta}\Im\left(f(x+i\delta
u)\right)u+\frac{n}{\delta}\xi u,\quad u\sim\mathbb{S}^{n-1}$ (3.6)

for some noise term $\xi\in\Xi$. In contrast to the noise-free setting in
[JYK21], equation (3.6) immediately reveals the delicacy in selecting
$\delta\in\mathbb{R}_{>0}$. Note, the term $n/\delta$ follows from our choice
to average over $\mathbb{B}^{n}$, i.e., by (3.1). Below we will clarify that
this term, and thereby the offset due to the noise, cannot be decreased by any
other choice of solid. In that sense, $\mathbb{B}^{n}$ is geometrically
optimal. We will use (3.6) in gradient descent algorithms of the form
$x_{k+1}=x_{k}-\mu_{k}g_{\delta_{k}}(x_{k})$, as detailed in Algorithm 4.1 (a)
and Algorithm 4.1 (b), for $\mu_{k}\in\mathbb{R}$ a stepsize and
$\delta_{k}\in\mathbb{R}_{>0}$ the smoothing parameter.

The next assumption on the (computational) noise will be assumed throughout.

######  Assumption 3.5 (Independence).

The random variable $\xi$ is drawn independently of $u\sim\mathbb{S}^{n-1}$.

######  Proposition 3.6 (Gradient approximation quality [JYK21, Proposition
3.4]).

Let $f\in C^{\omega,2}_{L_{2}(f)}(\mathcal{D})$ with $L_{2}(f)\geq 0$ satisfy
Assumption 3.1 for some $\bar{\delta}\in(0,1)$. Then, for any fixed
$x\in\mathcal{D}$ and $\kappa\in(0,1)$ there is a constant $C_{1}\geq 0$,
vanishing with $L_{2}(f)$, such that

$\displaystyle\|\nabla f_{\delta}(x)-\nabla f(x)\|_{2}$
$\displaystyle\leq{C_{1}n\delta^{2}}\quad\forall\delta\in(0,\kappa\bar{\delta}].$
(3.7)

We see that the simple single-point approach allows for an error of the form
$O(\delta^{2})$ which is what can be commonly achieved using central-
difference multi-points methods cf. [NS17].

From (3.7) it appears that (3.6) is potentially a biased gradient estimator.
Consider the special case of $f$ being quadratic (see Figure 2.2 for a
visualization). In that case, $\nabla f_{\delta}=\nabla f$, that is, the
estimator is unbiased: $\mathbb{E}[g_{\delta}]=\nabla f$. This property will
be exploited in Section 4.3.

In general, however, there will be a bias, controlled in part by selecting the
sequence $\\{\delta_{k}\\}_{k\geq 1}$ and unfortunately, a fixed bias
prohibits (local) convergence in general [AS21]. However, by looking at (3.6),
it can be shown that to overcome this, a selection of $\\{\mu_{k}\\}_{k\geq
1}$ and $\\{\delta_{k}\\}_{k\geq 1}$ should satisfy the following;

  1. (i)

As $\mu_{k}=\Theta(k^{-1})$ [RSS12], for fixed $\delta_{k}=\delta>0$ a bias
term prevails of the form $\sum^{K}_{k=1}\mu_{k}\delta=O(\log(K)+1)$. This can
be avoided by selecting $\delta_{k}$ to be asymptotically vanishing.

  2. (ii)

However, as the data is noisy, a term of the form $\mu_{k}/\delta_{k}$ also
accumulates. As such, by (i) $\delta_{k}\to 0$, but slower than $\mu_{k}\to
0$.

With this in mind we see that when $\mathbb{E}[g_{\delta}]\neq\nabla f$
zeroth-order optimization algorithms resort to selecting the smoothing-
parameter sequence $\\{\delta_{k}\\}_{k\geq 1}$ such that $\delta_{1}$
converges to $0$ sufficiently slow, cf. [NG21, Theorem 1], [BG21, Theorem 3].
See also [Fab71], [Spa05, Chapter 6], [WZS21, Assumption 1] for similar
assumptions from the stochastic approximation viewpoint. Motivated by the
observation that $\delta_{k}\to 0$ is necessary for an abundance of
algorithms, this work provides a framework that can handle this requirement
numerically. That means, a framework where $\delta_{k}$ can be made
arbitrarily small888Up to what the machine at hand can produce, usually
$2^{-1023}\approx 10^{-308}$..

At last, to characterize the effectiveness of our algorithms, we need to bound
the second moment of the estimator (3.6). We observe the same attractive
property as highlighted in [NS17], there is no need to assume boundedness of
the second moment of our stochastic estimator, cf. [RSS12]. As we allow for
computational noise, the bound will differ slightly from the result in
[JYK21].

######  Lemma 3.7 (Estimator second moment).

Let $f\in C^{\omega,2}_{{L_{2}}(f)}(\mathcal{D})$ satisfy Assumption 3.1 for
some $\bar{\delta}\in(0,1)$ and ${L_{2}}(f)\geq 0$. Then, for any fixed
$x\in\mathcal{D}$, $\kappa\in(0,1)$ and $g_{\delta}(x)$ as in (3.6) there are
constants $C_{a},C_{b}\geq 0$, vanishing with $L_{2}(f)$, such that for any
$\delta\in(0,\kappa\bar{\delta}]$ one has

$\displaystyle\mathbb{E}_{u\sim\mathbb{S}^{n-1}}\left[\|g_{\delta}(x)\|_{2}^{2}\right]$
$\displaystyle\leq{C_{a}n^{2}\delta^{4}}+C_{b}n^{2}\delta^{2}\|\nabla
f(x)\|_{2}+n\|\nabla f(x)\|_{2}^{2}+\tfrac{n^{2}}{\delta^{2}}\sigma_{\xi}.$
(3.8)

###### Proof.

First, observe from Algorithm 4.1, Assumption 1.2 and Assumption 3.5 that

$\mathbb{E}_{u\sim\mathbb{S}^{n-1}}\left[\|g_{\delta}(x)\|_{2}^{2}\right]=\tfrac{n^{2}}{\delta^{2}}\mathbb{E}_{u\sim\mathbb{S}^{n-1}}\left[\left(\Im\left(f(x+i\delta
u)\right)\right)^{2}\right]+\tfrac{n^{2}}{\delta^{2}}\mathbb{E}_{\xi}[\xi^{2}].$

Then, the claim follows directly by the same reasoning as in [JYK21, Corollary
3.5]. ∎

As with standard gradient-descent, the more isotropic the level sets of the
objective are, the better. The common way to enforce this is by means of
changing the underlying metric via the Hessian, i.e., Newton’s method. With
this in mind, averaging over some solid ellipsoid might appear more beneficial
than over the ball. In the spirit of [HL14] and [Hu+16, Proposition 3,Lemma 4]
we generalize Lemma 3.3 to more generic solids and show—perhaps
unsurprisingly— that spherical smoothing is optimal in the sense that it
minimizes the offset due to noise in (3.8).

To be in line with Assumption 3.1 we assume that this generic solid
$\mathsf{M}$ is a subset of $(-1,1)^{n}$.

######  Lemma 3.8 (The gradient of the complex-step function for generic
solids).

Let $\mathsf{M}\subset(-1,1)^{n}\subset\mathbb{R}^{n}$ be diffeomorphic to
$\mathbb{B}^{n}$. Let $f\in C^{\omega}(\mathcal{D})$ satisfy Assumption 3.1
for some $\bar{\delta}\in(0,1)$, then, $f_{\delta,\mathsf{M}}$ as in

$f_{\delta,\mathsf{M}}(x)=\mathbb{E}_{v\sim\mathsf{M}}\left[\Re(f(x+i\delta
v)\right]$ (3.9a) is differentiable and for any $x\in\mathcal{D}$ we have for
any $\delta\in(0,\bar{\delta})$ $\nabla
f_{\delta,\mathsf{M}}(x)=\frac{\mathrm{vol}(\delta\partial\mathsf{M})}{\mathrm{vol}(\delta\mathsf{M})}\cdot\mathbb{E}_{u\sim\partial\mathsf{M}}\left[\Im\left(f(x+i\delta
u)\right)N(u)\right].$ (3.9b)

for $N(u)$ a unit normal in $T^{\perp}_{u}\partial\mathsf{M}$.

###### Proof.

As $\mathsf{M}\subset\mathbb{R}^{n}$ is a compact oriented manifold with
boundary, we can appeal to the Divergence theorem [Lee13, Theorem 16.32]
(under the Euclidean metric), which states that for any smooth vector field
$X$ on $\mathsf{M}$ one has

$\int_{\mathsf{M}}\mathrm{div}(X(v))dV(v)=\int_{\partial\mathsf{M}}\langle
X(u),N(u)\rangle dV(u),$ (3.10)

for $N$ denoting the unit normal vector (field) along $\partial\mathsf{M}$.
That is, $\mathbb{R}^{n}=T_{p}\partial\mathsf{M}\oplus
T^{\perp}_{p}\partial\mathsf{M}$ for all $p\in\partial\mathsf{M}$ and $N(p)\in
T^{\perp}_{p}\partial\mathsf{M}$.

Using the same reasoning as for example in [JYK21], since one can select
$X=f\cdot C$ for $C$ some constant vector field on $\mathsf{M}$, then, as
$\mathrm{div}(C)=0$ and we can select $C$ to be aligned with any coordinate
axis, (3.10) implies that

$\int_{\mathsf{M}}\nabla f(v)dV(v)=\int_{\partial\mathsf{M}}f(u)N(u)dV(u).$
(3.11)

Note, $\nabla f(v)$ is well defined as $\mathsf{M}\subset\mathbb{R}^{n}$ is
diffeomorphic to $\mathbb{B}^{n}$.

Now we obtain the generalization of the result in [JYK21], that is, by
compactness, the Dominated Convergence theorem [Fol99, Section 2.3], the
Divergence theorem (3.10) and the Cauchy-Riemann equations [Kra00] we get

$\displaystyle\nabla_{x}\int_{\delta\mathsf{M}}\Re\left(f(x+iv)\right)\mathrm{d}V(v){=}\int_{\delta\partial\mathsf{M}}\Im\left(f(x+iu)\right)N(u)\mathrm{d}V(u),$

e.g., see [JYK21] for more on this line of reasoning. Then, due to the
distributional assumption (uniformity), we write

$\displaystyle f_{\delta,\mathsf{M}}(x)=$
$\displaystyle\mathbb{E}_{v\sim\mathsf{M}}\left[\Re\left(f(x+i\delta
v)\right)\right]=\frac{1}{\mathrm{vol}(\delta\mathsf{M})}\int_{\delta\mathsf{M}}\Re\left(f(x+iv)\right)\mathrm{d}V(v),$

and similarly,

$\displaystyle\mathbb{E}_{u\sim\partial\mathsf{M}}\left[\Im\left(f(x+i\delta
u)\right)N(u)\right]=\frac{1}{\mathrm{vol}(\delta\partial\mathsf{M})}\int_{\delta\partial\mathsf{M}}\Im\left(f(x+iu)\right)N(u)\mathrm{d}V(u).$

Combining it all yields (3.9b). ∎

As $N(u)\in T^{\perp}_{u}\partial\mathsf{M}$ is a unit vector, the offset term
in the variance (3.8) is minimized when we select $\mathsf{M}$ as

$\operatorname*{argmin}_{\mathsf{M}\in\mathscr{M}}\frac{\mathrm{vol}(\delta\partial\mathsf{M})}{\mathrm{vol}(\delta\mathsf{M})},$
(3.12)

where $\mathscr{M}$ is the set of manifolds diffeomorphic to $\mathbb{B}^{n}$
and $\delta\in\mathbb{R}_{>0}$. To retrieve the optimizer, consider the
isoperimetric inequality in $\mathbb{R}^{n}$ [Oss78] which implies that
$\mathsf{M}^{\star}=\mathbb{B}^{n}$ is optimal in the sense of (3.12).

To get (the complex-step version of) [HL14, Corollary 6] from Lemma 3.8, let
$\mathcal{E}^{n}_{Q}=\\{x\in\mathbb{R}^{n}:\langle Q^{-1}x,x\rangle\leq 1\\}$
for some $Q\in\mathcal{S}^{n}_{\succ 0}$. Now,
$T_{p}\partial\mathcal{E}^{n}_{Q}=\\{v\in\mathbb{R}^{n}:\langle
Q^{-1}p,v\rangle=0\\}$. As $\mathcal{E}^{n}_{Q}=Q^{1/2}\mathbb{B}^{n}$ one can
write

$f_{\delta,\mathcal{E}^{n}_{Q}}(x)=\mathbb{E}_{v\sim\mathcal{E}^{n}_{Q}}\left[f(x+i\delta
v)\right]=\mathbb{E}_{v\sim\mathbb{B}^{n}}\left[f(x+i\delta Q^{1/2}v)\right].$
(3.13a) Via the rightmost term in (3.13a) and the proof of Lemma 3.8 it
follows immediately that $\nabla
f_{\delta,\mathcal{E}^{n}_{Q}}(x)=\mathbb{E}_{u\sim\mathbb{S}^{n-1}}\frac{n}{\delta}\left[f(x+i\delta
Q^{1/2}u)Q^{-1/2}u\right].$ (3.13b)

Equivalently, one can directly appeal to (3.9b). However, here one needs to
appeal to the isoperimetric ratio for ellipsoids [Riv07].

At last, we provide further comments on applicability. The complex-step
derivative appears in a host of numerical applications, most notably, it is
reported in [CH04, Page 44] that a value of $\delta=10^{-100}$ is successfully
used in National Physical Laboratory software. In the context of zeroth-order
optimization, due to the complex-lifting, i.e., we need $f(x+i\delta u)$, we
cannot use immediately use physical measurement data, but we can work with any
simulation-based system or data that admits a complex representation. A few
areas of application are

  1. (i)

Simulation-based optimization, e.g., reinforcement learning and PDE-
constrained optimization, see also [WS21] and Example 5.3;

  2. (ii)

Privacy-sensitive optimization, e.g., the objective is known, but not to
everybody;

  3. (iii)

Black-box objective, e.g., $f(x)$ has been implemented in deprecated software,
see also [NS17].

##  4 Strongly convex imaginary zeroth-order optimization

In this section we will utilize the imaginary gradient estimator $g_{\delta}$
as given by (3.6) in the context of zeroth-order optimization algorithms. We
will not focus on fully generic convex optimization problems as the flat parts
of real-analytic convex functions must have measure zero [Kra00, JYK21].
Hence, without too much loss of generality we omit convex functions which are
not strongly convex999Future work will highlight the intimate relation between
convex and strongly convex functions under the assumption that both are real
analytic.. See also [KSST09] for more on strong-convexity in the context of
generalization.

In this section we relax some of the assumptions in [JYK21], not only can we
handle computational noise, the algorithms demand less knowledge of the
problem compared to other work. This is possible by introducing a time-varying
stepsize and a construction very much in line with [APT20]. In fact, recall
from [RSS12] that $\mu_{k}=\Theta(k^{-1})$ to allow for optimal rates. The
edge our results have, however, over these existing works is that our sequence
of smoothing parameters $\\{\delta_{k}\\}_{k\geq 1}$ is never catastrophic.

The generic algorithm for the unconstrained case is detailed in Algorithm 4.1
(a), i.e., for $\mathcal{X}=\mathcal{D}$.

1: Input: initial iterate $x_{1}\in\mathcal{X}$, stepsizes
$\\{\mu_{k}\\}_{k\geq 1}$, smoothing parameters $\\{\delta_{k}\\}_{k\geq 1}$.

2: for $k=1,2,\ldots,K$ do

3: generate random $u_{k}\sim\mathbb{S}^{n-1}$

4: obtain noisy estimate
$g_{\delta_{k}}(x_{k})=\frac{n}{\delta_{k}}\Im\left(f(x_{k}+i\delta_{k}u_{k})\right)u_{k}+\frac{n}{\delta_{k}}\xi_{k}u_{k}$

5: set $x_{k+1}=\Pi_{\mathcal{X}}\left(x_{k}-\mu_{k}\cdot
g_{\delta_{k}}(x_{k})\right)$

6: end for

Algorithm 4.1 Imaginary zeroth-order optimization:
(a) unconstrained $\mathcal{X}=\mathcal{D}$ and (b) constrained
$\mathcal{X}=\mathcal{K}\subset\mathcal{D}$.

Given a compact (possibly non-convex) set $\mathcal{K}\subset\mathcal{D}$, the
algorithm for the constrained case is detailed in Algorithm 4.1 (b) i.e., for
$\mathcal{X}=\mathcal{K}$. Here,
$\Pi_{\mathcal{K}}:\mathcal{D}\subseteq\mathbb{R}^{n}\to\mathcal{K}$ denotes
the projection operator.

Note, in our algorithms we will assume that we can sample in a small
$\delta$-neighbourhood contained in $\mathcal{D}\setminus\mathcal{K}$. As a
key application of the proposed framework is simulation-based optimization
this is deemed justifiable. Having access to a projection operator
$\Pi_{\mathcal{K}}$, we will assume nothing more than feasibility regarding
the initial condition $x_{1}$.

###  4.1 Strong convexity

In this part we consider the setting of $f\in C^{\omega}(\mathcal{D})$ being
$\tau(f)$-strongly convex over $\mathcal{D}$, i.e., there is some $\tau(f)>0$
such that

$f(y)\geq f(x)+\langle\nabla
f(x),y-x\rangle+\tfrac{1}{2}\tau(f)\|y-x\|_{2}^{2},\quad\forall
x,y\in\mathcal{D}.$ (4.1)

In particular (4.1) implies that for $\mathcal{D}$ such that
$x^{\star}\in\mathrm{int}(\mathcal{D})$ one has

$f(x)-f(x^{\star})\geq\tfrac{1}{2}\tau(f)\|x-x^{\star}\|_{2}^{2},\quad\forall
x\in\mathcal{D}.$ (4.2)

If additionally $f\in C^{\omega,1}_{L_{1}(f)}$, then by $\|\nabla
f(x)\|_{2}^{2}\geq 2\tau(f)(f(x)-f(x^{\star}))$ one has

$\tau(f)\|x-x^{\star}\|_{2}\leq\|\nabla f(x)\|_{2}\leq
L_{1}(f)\|x-x^{\star}\|_{2}.$ (4.3)

In contrast to [JYK21], our algorithms “only” demand knowledge of the strong-
convexity parameter. In Section 4.5 we mention how one could estimate
$\tau(f)$.

###  4.2 Generic convergence rates

As in [APT20], we start with the constrained case.

######  Theorem 4.2 (Convergence rate of Algorithm 4.1 (b) with noise).

Let $f\in C^{\omega}(\mathcal{D})$ be a $\tau(f)$-strongly convex function
satisfying Assumption 3.1 for some $\bar{\delta}\in(0,1)$ and let
$\mathcal{K}\subset\mathcal{D}$ be a compact convex set. Suppose that $f$ has
a Lipschitz Hessian over $\mathcal{K}$, that is, (1.3) holds for a non-zero
constant $L_{2}(f)$. Let $\\{x_{k}\\}_{k\geq 1}$ be the sequence of iterates
generated by Algorithm 4.1 (b) with stepsize $\mu_{k}={2}/({\tau(f)k})$ and
the sequence of smoothing parameters defined for all $k\geq 1$ by
$\delta_{k}=\delta k^{-1/6}$ with $\delta\in(0,\kappa\bar{\delta}]$ for some
$\kappa\in(0,1)$. Then, if the oracle satisfies Assumption 1.2, the uniformly-
averaged iterate $\bar{x}_{K}=K^{-1}\sum^{K}_{k=1}x_{k}$ achieves the
optimization error

$\mathbb{E}[f(\bar{x}_{K})-f(x^{\star})]\leq\widetilde{O}\left(\frac{n^{2}}{\tau(f)}\delta^{-\tfrac{1}{3}}\sigma_{\xi}K^{-\tfrac{2}{3}}\right).$

###### Proof.

We mainly follow [APT20]. To that end, let $\sup_{x\in\mathcal{K}}\|\nabla
f(x)\|_{2}\leq G$. As $\mathcal{K}$ is convex and compact we have by the
properties of the operator $\Pi_{\mathcal{K}}$ that
$\|x_{k+1}-x^{\star}\|_{2}^{2}\leq\|x_{k}-\mu_{k}g_{\delta_{k}}(x_{k})-x^{\star}\|_{2}^{2}$.
This can be written as conveniently as

$\langle
g_{\delta_{k}}(x_{k}),x_{k}-x^{\star}\rangle\leq\tfrac{1}{2\mu_{k}}\left(\|x_{k}-x^{\star}\|_{2}^{2}-\|x_{k+1}-x^{\star}\|_{2}^{2}\right)+\tfrac{\mu_{k}}{2}\|g_{\delta_{k}}(x_{k})\|_{2}^{2}.$
(4.4)

After reordering the standard strong $\tau(f)$-convexity expression, one
obtains

$f(x_{k})-f(x^{\star})\leq\langle\nabla
f(x_{k}),x_{k}-x^{\star}\rangle-\tfrac{\tau(f)}{2}\|x_{k}-x^{\star}\|_{2}^{2}.$
(4.5)

Set $a_{k}=\|x_{k}-x^{\star}\|_{2}^{2}$, then, an application of the Cauchy-
Schwarz inequality after combining (4.4) with (4.5) and taking the expectation
over $u_{k}$ and $\xi_{k}$ conditioned on $x_{k}$ yields

$\displaystyle\begin{array}[]{lcl}\mathbb{E}[f(x_{k})-f(x^{\star})|x_{k}]&\leq&\left\|\mathbb{E}[g_{\delta_{k}}(x_{k})|x_{k}]-\nabla
f(x_{k})\right\|_{2}\|x_{k}-x^{\star}\|_{2}+\tfrac{1}{2\mu_{k}}\mathbb{E}[a_{k}-a_{k+1}|x_{k}]\\\
&&\tfrac{\mu_{k}}{2}\mathbb{E}[\|g_{\delta_{k}}(x_{k})\|_{2}^{2}|x_{k}]-\tfrac{\tau(f)}{2}\mathbb{E}[a_{k}|x_{k}]\\\
&\overset{\eqref{equ:grad:delta:approx:error}}{\leq}&C_{1}n\delta_{k}^{2}\|x_{k}-x^{\star}\|_{2}+\tfrac{1}{2\mu_{k}}\mathbb{E}[a_{k}-a_{k+1}|x_{k}]\\\
&&\tfrac{\mu_{k}}{2}\mathbb{E}[\|g_{\delta_{k}}(x_{k})\|_{2}^{2}|x_{k}]-\tfrac{\tau(f)}{2}\mathbb{E}[a_{k}|x_{k}],\end{array}$

for some $C_{1}>0$. Now, use $ab\leq\tfrac{1}{2}(a^{2}+b^{2})$, in particular
$ab\leq\tfrac{1}{2}(\gamma a^{2}+\gamma^{-1}b^{2})$ for $\gamma\neq 0$, to
construct

$n\delta_{k}^{2}\|x_{k}-x^{\star}\|_{2}\leq\tfrac{1}{2}\left(\tfrac{2C_{1}}{\tau(f)}n^{2}\delta_{k}^{4}+\tfrac{\tau(f)}{2C_{1}}\|x_{k}-x^{\star}\|_{2}^{2}\right).$

Next, take the expectation over $x_{k}$ and let $r_{k}=\mathbb{E}[a_{k}]$ such
that we can write

$\mathbb{E}[f(x_{k})-f(x^{\star})]\leq\tfrac{1}{2\mu_{k}}(r_{k}-r_{k+1})-\tfrac{\tau(f)}{4}r_{k}+\tfrac{1}{\tau(f)}C_{1}^{2}n^{2}\delta_{k}^{4}+\tfrac{\mu_{k}}{2}\mathbb{E}[\|g_{\delta_{k}}(x_{k})\|_{2}^{2}].$
(4.6)

Summing (4.6) over $k$ yields

$\displaystyle\textstyle\sum^{K}_{k=1}\mathbb{E}[f(x_{k})-f(x^{\star})]\leq\tfrac{1}{2}$
$\displaystyle\textstyle\sum^{K}_{k=1}\left(\tfrac{1}{\mu_{k}}(r_{k}-r_{k+1})-\tfrac{\tau(f)}{2}r_{k}\right)$
$\displaystyle+$
$\displaystyle\textstyle\sum^{K}_{k=1}\left(\tfrac{1}{\tau(f)}C_{1}^{2}n^{2}\delta_{k}^{4}+\tfrac{\mu_{k}}{2}\mathbb{E}[\|g_{\delta_{k}}(x_{k})\|_{2}^{2}]\right).$

As we selected $\mu_{k}=2/(\tau(f)k)$ we can simplify the above by using the
same reasoning as in [APT20], that is

$\displaystyle\textstyle\sum^{K}_{k=1}\left(\tfrac{1}{\mu_{k}}(r_{k}-r_{k+1})-\tfrac{\tau(f)}{2}r_{k}\right)\leq
r_{1}\left(\tfrac{1}{\mu_{1}}-\tfrac{\tau(f)}{2}\right)+\textstyle\sum^{K}_{k=2}r_{k}\left(\tfrac{1}{\mu_{k}}-\tfrac{1}{\mu_{k-1}}-\tfrac{\tau(f)}{2}\right)=0.$

Note that we rely on the $\tau(f)$-strong convexity. Using the observation
from above and plugging in the stepsize $\mu_{k}$ elsewhere yields by (3.8)

$\displaystyle\textstyle\sum^{K}_{k=1}\mathbb{E}[f(x_{k})-f(x^{\star})]$
$\displaystyle\leq\tfrac{1}{\tau(f)}\textstyle\sum^{K}_{k=1}\left(C_{1}^{2}n^{2}\delta_{k}^{4}+\tfrac{1}{k}\mathbb{E}[\|g_{\delta_{k}}(x_{k})\|_{2}^{2}]\right)$
$\displaystyle\leq\tfrac{n^{2}}{\tau(f)}\textstyle\sum^{K}_{k=1}\left(C_{1}^{2}\delta_{k}^{4}+\tfrac{1}{k}\left[{C_{2}\delta_{k}^{4}}+C_{3}\delta_{k}^{2}\|\nabla
f(x_{k})\|_{2}+\tfrac{1}{n}\|\nabla
f(x_{k})\|_{2}^{2}+\frac{1}{\delta_{k}^{2}}\sigma_{\xi}\right]\right),$

for some $C_{2},C_{3}>0$. Now, minimizing over $\\{\delta_{k}\\}_{k}$ is
possible but yields smoothing parameters as a function of unknown constants.
Instead, we retain the “optimal” root101010Let $a,b\in\mathbb{R}_{>0}$, then,
see that
$(b/(2a))^{\tfrac{1}{6}}=\operatorname*{argmin}_{\delta\in\mathbb{R}_{\geq
0}}\\{a\delta^{4}+b\tfrac{1}{\delta^{2}}\\}$. and propose

$\widetilde{\delta}_{k}=\left(\frac{\alpha\sigma_{\xi}}{k}\right)^{\tfrac{1}{6}},$

for some $\alpha\in(0,1)$ to be specified. Using this smoothing parameter
sequence, that is, $\delta_{k}=\widetilde{\delta}_{k}$, together with
$\sum^{K}_{k=1}k^{-1}\leq 1+\log(K)$ (Lemma A.1) yields

$\displaystyle\begin{array}[]{lcl}\textstyle\sum^{K}_{k=1}\mathbb{E}[f(x_{k})-f(x^{\star})]&\leq&\tfrac{n^{2}}{\tau(f)}\textstyle\sum^{K}_{k=1}\left(C_{1}^{2}\left(\frac{\alpha\sigma_{\xi}}{k}\right)^{\tfrac{2}{3}}+\tfrac{1}{k}\left[{C_{2}\left(\frac{\alpha\sigma_{\xi}}{k}\right)^{\tfrac{2}{3}}}+\left(\frac{\alpha\sigma_{\xi}}{k}\right)^{-\tfrac{1}{3}}\sigma_{\xi}\right]\right)\\\
&&+\tfrac{n}{\tau(f)}G^{2}(1+\log(K))+\textstyle\tfrac{n^{2}}{\tau(f)}GC_{3}\sum^{K}_{k=1}\tfrac{1}{k}\left(\frac{\alpha\sigma_{\xi}}{k}\right)^{\tfrac{1}{3}}\\\
&=&\tfrac{n^{2}}{\tau(f)}\textstyle\sum^{K}_{k=1}\left(C_{1}^{2}\left(\frac{\alpha\sigma_{\xi}}{k}\right)^{\tfrac{2}{3}}+\tfrac{1}{k}\left[{C_{2}\left(\frac{\alpha\sigma_{\xi}}{k}\right)^{\tfrac{2}{3}}}+k^{\tfrac{1}{3}}\sigma_{\xi}^{\tfrac{2}{3}}\alpha^{-\tfrac{1}{3}}\right]\right)\\\
&&+\tfrac{n}{\tau(f)}G^{2}(1+\log(K))+\textstyle\tfrac{n^{2}}{\tau(f)}GC_{3}\sum^{K}_{k=1}k^{-\tfrac{2}{3}}(\alpha\sigma_{\xi})^{\tfrac{1}{3}}\\\
&\leq&\tfrac{n^{2}}{\tau(f)}\textstyle\sum^{K}_{k=1}C_{4}k^{-\tfrac{2}{3}}\sigma_{\xi}^{\tfrac{2}{3}}\alpha^{-\tfrac{1}{3}}\\\
&&+\tfrac{n}{\tau(f)}G^{2}(1+\log(K))+\textstyle\tfrac{n^{2}}{\tau(f)}GC_{3}\sum^{K}_{k=1}k^{-\tfrac{2}{3}}(\alpha\sigma_{\xi})^{\tfrac{1}{3}}.\end{array}$

Now, as $\sum^{K}_{k=1}k^{-\tfrac{2}{3}}\leq 3K^{\tfrac{1}{3}}$ (Lemma A.2) we
can continue and write

$\displaystyle\begin{array}[]{lcl}\textstyle\sum^{K}_{k=1}\mathbb{E}[f(x_{k})-f(x^{\star})]&\leq&\tfrac{n^{2}}{\tau(f)}C_{5}K^{\tfrac{1}{3}}\sigma_{\xi}^{\tfrac{2}{3}}\alpha^{-\tfrac{1}{3}}+\tfrac{n}{\tau(f)}G^{2}(1+\log(K))+\textstyle\tfrac{n^{2}}{\tau(f)}GC_{6}K^{\tfrac{1}{3}}(\alpha\sigma_{\xi})^{\tfrac{1}{3}}.\end{array}$

and as such we obtain the optimization error

$\displaystyle\begin{array}[]{lcl}\mathbb{E}[f(\bar{x}_{K})-f(x^{\star})]&\leq&\tfrac{n^{2}}{\tau(f)}C_{5}K^{-\tfrac{2}{3}}\sigma_{\xi}^{\tfrac{2}{3}}\alpha^{-\tfrac{1}{3}}+\tfrac{n}{\tau(f)}G^{2}K^{-1}(1+\log(K))+\textstyle\tfrac{n^{2}}{\tau(f)}GC_{6}K^{-\tfrac{2}{3}}(\alpha\sigma_{\xi})^{\tfrac{1}{3}}.\end{array}$

As $\alpha\in(0,1)$ was arbitrary, we can set $\delta=\alpha\sigma_{\xi}$ such
that $\delta_{k}=\delta k^{-\tfrac{1}{6}}$ for some
$\delta\in(0,\bar{\delta})$. ∎

The edge Theorem 4.2 has over existing work is that the requested sequence
$\\{\delta_{k}\\}_{k\geq 1}$ can always be safely implemented. With respect to
optimality, we highlight a general method to pass from $\widetilde{O}(\cdot)$
to $O(\cdot)$ complexities.

######  Remark 4.3 (Removing the logarithmic term).

One can appeal to $\alpha$-suffix averaging as proposed in [RSS12] to remove
the logarithmic term. This is achieved by averaged estimates of the form
$\widetilde{x}_{K}=\tfrac{2}{K}\sum^{K}_{k=K/2+1}x_{k}$ and follows from
$\sum^{T}_{t=(1-\alpha)T+1}\tfrac{1}{t}\leq\log(1/(1-\alpha))$ for
$\alpha\in(0,1)$ such that $\alpha T,(1-\alpha)T\in\mathbb{Z}$. As the
implementation of $\widetilde{x}_{K}$ is not always easier or more efficient
than $\bar{x}_{K}$, the uniformly-averaged estimator remains competitive
despite the slower rate.

Next we consider the unconstrained case. Here, we cannot appeal to an uniform
bound on $\nabla f(x)$. Instead, we use the idea from [APT20, Theorem 3.2] and
bound a subset of iterates before strong-convexity kicks in. In practise, when
$\tau(f)$ is small, the first few stepsizes will be relatively large and can
lead to overflow. In some sense one could interpret this as some restarting
mechanism.

######  Theorem 4.4 (Convergence rate of Algorithm 4.1 (a) with noise).

Let $f\in C^{\omega}(\mathcal{D})$ be a $\tau(f)$-strongly convex function
satisfying Assumption 3.1 for some $\bar{\delta}\in(0,1)$ with
$x^{\star}\in\mathrm{int}(\mathcal{D})$. Suppose that $f$ has a Lipschitz
gradient and Hessian, that is, (1.2) and (1.3) hold, for non-zero constants
$L_{1}(f)$ and $L_{2}(f)$, respectively. Let $\\{x_{k}\\}_{k\geq 1}$ be the
sequence of iterates generated by Algorithm 4.1 (a) for

$\displaystyle\begin{array}[]{llll}&\mu_{k}=\tfrac{1}{\tau(f)K},&\delta_{k}=\delta
K^{-\tfrac{1}{6}},&k=1,\dots,K_{0},\\\
&\mu_{k}=\tfrac{2}{\tau(f)k},&\delta_{k}=\delta
k^{-\tfrac{1}{6}},&k=K_{0}+1,\dots,K,\end{array}$

with $K_{0}=\left\lfloor\tfrac{8n^{2}L_{1}(f)^{2}}{\tau(f)^{2}}\right\rfloor$
and $\delta\in(0,\kappa\bar{\delta}]$ for some $\kappa\in(0,1)$. Then, if the
oracle satisfies Assumption 1.2 and $K\geq 2K_{0}$ we incur for
$\bar{x}_{K_{0},K}=\tfrac{1}{K-K_{0}}\sum^{K}_{k=K_{0}+1}x_{k}$ the
optimization error

$\displaystyle\mathbb{E}[f(\bar{x}_{K_{0},K})-f(x^{\star})]$
$\displaystyle\leq
O\left(\frac{n^{2}L_{1}(f)^{2}}{\tau(f)}\|x_{1}-x^{\star}\|_{2}^{2}K^{-1}\right)+O\left(\frac{n^{2}\sigma_{\xi}}{\tau(f)\delta^{2}}K^{-\tfrac{2}{3}}\right).$
(4.7)

###### Proof.

The proof will be similar to that of [APT20, Theorem 3.2]. Again, set
$a_{k}=\|x_{k}-x^{\star}\|_{2}^{2}$, then, as in the proof of Theorem 4.2

$\displaystyle\begin{array}[]{lcl}\mathbb{E}[f(x_{k})-f(x^{\star})|x_{k}]&\leq&\left\|\mathbb{E}[g_{\delta_{k}}(x_{k})|x_{k}]-\nabla
f(x_{k})\right\|_{2}\|x_{k}-x^{\star}\|_{2}+\tfrac{1}{2\mu_{k}}\mathbb{E}[a_{k}-a_{k+1}|x_{k}]\\\
&&\tfrac{\mu_{k}}{2}\mathbb{E}[\|g_{\delta_{k}}(x_{k})\|_{2}^{2}|x_{k}]-\tfrac{\tau(f)}{2}\mathbb{E}[a_{k}|x_{k}]\\\
&\overset{\eqref{equ:grad:delta:approx:error}}{\leq}&C_{1}n\delta_{k}^{2}\|x_{k}-x^{\star}\|_{2}+\tfrac{1}{2\mu_{k}}\mathbb{E}[a_{k}-a_{k+1}|x_{k}]\\\
&&\tfrac{\mu_{k}}{2}\mathbb{E}[\|g_{\delta_{k}}(x_{k})\|_{2}^{2}|x_{k}]-\tfrac{\tau(f)}{2}\mathbb{E}[a_{k}|x_{k}].\end{array}$

Now, use $ab\leq\tfrac{1}{2}(a^{2}+b^{2})$ together with $\tau(f)$-strong
convexity, i.e., (4.2), to construct

$n\delta_{k}^{2}\|x_{k}-x^{\star}\|_{2}\leq\tfrac{1}{2}\left(\tfrac{2C_{1}}{\tau(f)}n^{2}\delta_{k}^{4}+\tfrac{\tau(f)}{2C_{1}}\|x_{k}-x^{\star}\|_{2}^{2}\right)\leq\tfrac{C_{1}}{\tau(f)}n^{2}\delta_{k}^{4}+\tfrac{1}{2C_{1}}(f(x_{k})-f(x^{\star})).$

Next, let $r_{k}=\mathbb{E}[a_{k}]$ such that by $\|\nabla
f(x_{k})\|_{2}^{2}\leq L_{1}(f)^{2}\|x_{k}-x^{\star}\|_{2}^{2}$ we can write

$\displaystyle\begin{array}[]{lcl}\mathbb{E}[f(x_{k})-f(x^{\star})]&\leq&\tfrac{1}{\mu_{k}}(r_{k}-r_{k+1})-\tau(f)r_{k}+\tfrac{2}{\tau(f)}C_{1}^{2}n^{2}\delta_{k}^{4}+{\mu_{k}}\mathbb{E}[\|g_{\delta_{k}}(x_{k})\|_{2}^{2}]\\\
&\overset{\eqref{equ:g:var}}{\leq}&\tfrac{1}{\mu_{k}}(r_{k}-r_{k+1})-{\tau(f)}r_{k}+\tfrac{2}{\tau(f)}C_{1}^{2}n^{2}\delta_{k}^{4}\\\
&&+{\mu_{k}}\left({C_{2}n^{2}\delta_{k}^{4}}+C_{3}n^{2}\delta_{k}^{4}+2n^{2}L_{1}(f)^{2}r_{k}+\tfrac{n^{2}}{\delta_{k}^{2}}\sigma_{\xi}\right)\end{array}$
(4.8)

Where in the last step we used

$\delta_{k}^{2}\|x_{k}-x^{\star}\|_{2}\leq\tfrac{1}{2}\left(\tfrac{C_{3}}{2L_{1}(f)}\delta_{k}^{4}+\tfrac{2L_{1}(f)}{C_{3}}\|x_{k}-x^{\star}\|_{2}^{2}\right)$

to rewrite (3.8).

Now we use the step- and smoothingsize for $k=1,\dots,K_{0}$, that is,
$\mu_{k}=1/(\tau(f)K)$, $\delta_{k}=\delta K^{-\tfrac{1}{6}}$, and observe
that

$\displaystyle r_{k+1}$ $\displaystyle\leq
r_{k}-\tau(f)\mu_{k}r_{k}+\tfrac{2\mu_{k}}{\tau(f)}C_{1}^{2}n^{2}\delta_{k}^{4}+\mu_{k}^{2}\left({C_{2}n^{2}\delta_{k}^{4}}+C_{3}n^{2}\delta_{k}^{4}+2n^{2}L_{1}(f)^{2}r_{k}+\tfrac{n^{2}}{\delta_{k}^{2}}\sigma_{\xi}\right)$
$\displaystyle=\left(1-\tfrac{1}{K}+\tfrac{2L_{1}(f)^{2}}{(\tau(f)K)^{2}}n^{2}\right)r_{k}+\nu_{K}=a_{K}r_{k}+\nu_{K}$

for $a_{k}$ as between brackets and $\nu_{K}$ defined as

$\displaystyle\nu_{K}=$
$\displaystyle\tfrac{2}{\tau(f)^{2}K}C_{1}^{2}n^{2}\delta_{k}^{4}+\tfrac{1}{(\tau(f)K)^{2}}\left({C_{2}n^{2}\delta_{k}^{4}}+C_{3}n^{2}\delta_{k}^{4}+\tfrac{n^{2}}{\delta_{k}^{2}}\sigma_{\xi}\right)$
$\displaystyle\leq$
$\displaystyle\tfrac{1}{\tau(f)^{2}K}\left(n^{2}\delta^{4}C_{4}+\tfrac{n^{2}}{\delta^{2}}\sigma_{\xi}\right)K^{-\tfrac{2}{3}}.$

We now proceed with bounding $r_{K_{0}+1}$. As in [APT20], set

$q_{K}=1+\tfrac{2L_{1}(f)^{2}}{(\tau(f)K)^{2}}n^{2}$

by iterating over $r_{k}$ it follows from a geometric series argument that

$\displaystyle r_{K_{0}+1}\leq
a_{K}^{K_{0}}r_{1}+\textstyle\sum^{K_{0}-1}_{i=0}a_{K}^{i}\nu_{K}\leq\left(r_{1}+\tfrac{(\tau(f)K)^{2}}{2L_{1}(f)^{2}n^{2}}\nu_{K}\right)q_{K}^{K_{0}}.$

Now for $\left\lfloor\cdot\right\rfloor$ being the floor function, let $K_{0}$
be as in the theorem. Then, as $\log(1+x)\leq x$, on $\mathbb{R}_{\geq 0}$,
one has

$\displaystyle q_{K}^{K_{0}}=$
$\displaystyle\mathrm{exp}\left(K_{0}\log\left(1+\tfrac{2L_{1}(f)^{2}}{(\tau(f)K)^{2}}n^{2}\right)\right)\leq\mathrm{exp}\left(\tfrac{8n^{2}L_{1}(f)^{2}}{\tau(f)^{2}}\log\left(1+\tfrac{2L_{1}(f)^{2}}{(\tau(f)K)^{2}}n^{2}\right)\right)\leq\mathrm{exp}\left(\tfrac{16n^{4}L_{1}(f)^{4}}{\tau(f)^{4}K^{2}}\right).$

Fix any $\bar{C}\in(0,\tfrac{1}{32})$, when

$K=\sqrt{\tfrac{8n^{4}L_{1}(f)^{4}}{\tau(f)^{4}\bar{C}}}$

then $K\geq 2K_{0}$ and $q^{K_{0}}_{K}\leq e^{\bar{C}}=C_{5}$. As such,

$\displaystyle r_{K_{0}+1}\leq$
$\displaystyle\left(r_{1}+\tfrac{(\tau(f)K)^{2}}{2L_{1}(f)^{2}n^{2}}\nu_{K}\right)C_{5}$
$\displaystyle\leq$
$\displaystyle\left(r_{1}+\tfrac{(\tau(f)K)^{2}}{L_{1}(f)^{2}n^{2}}\tfrac{1}{\tau(f)^{2}K}\left(n^{2}\delta^{4}C_{4}+\tfrac{n^{2}}{\delta^{2}}\sigma_{\xi}\right)K^{-\tfrac{2}{3}}\right)C_{5}$
$\displaystyle=$
$\displaystyle\left(r_{1}+\tfrac{1}{L_{1}(f)^{2}n^{2}}\left(n^{2}\delta^{4}C_{4}+\tfrac{n^{2}}{\delta^{2}}\sigma_{\xi}\right)K^{\tfrac{1}{3}}\right)C_{5}.$

Now we return to our normal step- and smoothingsizes, that is
$\mu_{k}=2/(\tau(f)k)$, $\delta_{k}=\delta k^{-\tfrac{1}{6}}$, for $k\geq
K_{0}+1$. By plugging this into (4.8) we get

$\displaystyle(K-K_{0})\mathbb{E}[f(\bar{x}_{K_{0},K})-f(x^{\star})]$
$\displaystyle\leq\textstyle\sum^{K}_{k=K_{0}+1}\tfrac{\tau(f)k}{2}(r_{k}-r_{k+1})-{\tau(f)}r_{k}+\tfrac{4}{\tau(f)k}n^{2}L_{1}(f)^{2}r_{k}$
$\displaystyle+\textstyle\sum^{K}_{k=K_{0}+1}\tfrac{2}{\tau(f)k}C_{6}n^{2}\delta^{4}k^{-\tfrac{2}{3}}$
$\displaystyle+\textstyle\sum^{K}_{k=K_{0}+1}\tfrac{2}{\tau(f)k}\tfrac{n^{2}}{\delta^{2}}k^{\tfrac{1}{3}}\sigma_{\xi}$

By construction of $K_{0}$ we have that for $k\geq K_{0}+1$,
$\tau(f)/2\geq(4n^{2}L_{1}(f)^{2})/(\tau(f)k)$. Hence

$\displaystyle(K-K_{0})\mathbb{E}[f(\bar{x}_{K_{0},K})-f(x^{\star})]$
$\displaystyle\leq\tfrac{\tau(f)}{2}\left(\textstyle\sum^{K}_{k=K_{0}+1}k(r_{k}-r_{k+1})-r_{k}\right)+U_{K_{0},K}.$

where by Lemma A.2

$\displaystyle U_{K_{0},K}=$
$\displaystyle\textstyle\sum^{K}_{k=K_{0}+1}\tfrac{2}{\tau(f)k}C_{6}n^{2}\delta^{4}k^{-\tfrac{2}{3}}+\tfrac{2}{\tau(f)k}\tfrac{n^{2}}{\delta^{2}}k^{\tfrac{1}{3}}\sigma_{\xi}$
$\displaystyle\leq$
$\displaystyle\tfrac{2n^{2}}{\tau(f)}\left(C_{6}\delta^{4}+\tfrac{1}{\delta^{2}}\sigma_{\xi}\right)\textstyle\sum^{K}_{k=K_{0}+1}k^{-\tfrac{2}{3}}$
$\displaystyle\leq$
$\displaystyle\tfrac{2n^{2}}{\tau(f)}\left(C_{6}\delta^{4}+\tfrac{1}{\delta^{2}}\sigma_{\xi}\right)3K^{\tfrac{1}{3}}.$

As demonstrated in [APT20] (below), one can now construct the bound
$\textstyle\sum^{K}_{k=K_{0}+1}k(r_{k}-r_{k+1})-r_{k}\leq K_{0}r_{K_{0}+1}$
where the last term is exactly the term we could bound before. In combination
with the bound on $K_{0}$ itself, we find that

$\displaystyle(K-K_{0})\mathbb{E}[f(\bar{x}_{K_{0},K})-f(x^{\star})]$
$\displaystyle\leq\tfrac{\tau(f)}{2}\tfrac{8n^{2}L_{1}(f)^{2}}{\tau(f)^{2}}\left(r_{1}+\tfrac{1}{L_{1}(f)^{2}n^{2}}\left(n^{2}\delta^{4}C_{4}+\tfrac{n^{2}}{\delta^{2}}\sigma_{\xi}\right)K^{\tfrac{1}{3}}\right)C_{5}$
$\displaystyle+\tfrac{2n^{2}}{\tau(f)}\left(C_{6}\delta^{4}+\tfrac{1}{\delta^{2}}\sigma_{\xi}\right)3K^{\tfrac{1}{3}}.$

By our selection of $\bar{C}$ we have that $K\geq 2K_{0}$ and as such

$\displaystyle\mathbb{E}[f(\bar{x}_{K_{0},K})-f(x^{\star})]$
$\displaystyle\leq\tfrac{8n^{2}L_{1}(f)^{2}}{\tau(f)K}\left(r_{1}+\tfrac{1}{L_{1}(f)^{2}}\left(\delta^{4}C_{3}+\tfrac{1}{\delta^{2}}\sigma_{\xi}\right)K^{\tfrac{1}{3}}\right)C_{5}$
$\displaystyle+\tfrac{4n^{2}}{\tau(f)K}\left(C_{6}\delta^{4}+\tfrac{1}{\delta^{2}}\sigma_{\xi}\right)3K^{\tfrac{1}{3}}.$

Now, reordering terms yields (4.7). ∎

######  Remark 4.5 (On unconstrained anytime algorithms).

The unconstrained algorithms (Theorem 4.4 and Theorem 4.9) require the user to
pre-define the full length of the algorithm as the stepsize depends on $K$.
One can mitigate this by shifting the dependence on $k$, e.g., by using
$\mu_{k}=1/\big{(}\tau(f)(k+2K_{0})\big{)}$. Although the rate (4.7) remains
unaffected, this does come at the cost of potentially sacrificing progress in
the first $K_{0}$ steps of the algorithm. A detailled study is left for future
work.

###  4.3 Optimal convergence rates

Now we consider the special case of $f$ being quadratic. Here we improve upon
the previous section due to exploitation of the quadratic nature of $f$, that
is, by using $\Im(f(x+i\delta u))=\delta\langle\nabla f(x),u\rangle$ for any
$\delta>0$.

Better yet, we see that for quadratic functions we incur optimal regret.
Optimality can be shown along the lines of [Sha13], or along the lines of
[Aga+09] after observing that in the quadratic case the gradient estimator
$g_{\delta}(x)$ becomes an unbiased estimator for $\nabla f(x)$. The test
function used in [APT20] is smooth but unfortunately not analytic111111Section
6.2 highlights that this might not be an obstruction.. We start by providing
the bound from below.

######  Theorem 4.6 (Bound from below).

Any possibly randomized zeroth-order algorithm of fixed length $K$, applying
the estimator (3.6) under Assumption 1.2, cannot achieve a rate faster than

$\Omega\left(\frac{n^{2}}{\tau(f)K}\right),$

uniformly over all $\tau(f)$-strongly convex quadratic (real-analytic)
functions.

###### Proof.

We largely follow [Sha13, Theorem 3], but for the sake of completeness we
highlight the main arguments.

Recall that based on $x_{1},x_{2},\dots,x_{K}$, in particular the function
evaluations at those points, we compute some point ${x}^{\prime}_{K}$ (this
could be a non-uniform average estimator). In our case the function queries
correspond to $v_{k}=\Im(f(x_{k}+i\delta u))$ for some choice of $\delta>0$,
$u\in\mathbb{S}^{n-1}$ and with the possibility of being corrupted by additive
noise $\xi$.

Now, consider the $C^{\omega}$ function over $\mathbb{R}^{n}$

$f_{z}(x)=\tfrac{\tau}{2}\|x\|_{2}^{2}-\langle z,x\rangle.$ (4.9)

The unique minimizer of $f_{z}(x)$ is given by $x^{\star}=\tfrac{1}{\tau}z$.
Moreover, assume $z$ is drawn uniformly from $\\{-\nu,\nu\\}^{n}$ for some
$\nu$ that will be specified later. It follows from the strong
$\tau$-convexity of (4.9) that
$f_{z}(x)-f_{z}(x^{\star})\geq\tfrac{\tau}{2}\|x-\tfrac{1}{\tau}z\|_{2}^{2}$.
As such, for any randomized strategy

$\mathbb{E}[f_{z}(x^{\prime}_{K})-f_{z}(x^{\star})]\geq\tfrac{\tau}{2}\mathbb{E}[\|x^{\prime}_{K}-\tfrac{1}{\tau}z\|_{2}^{2}]=\tfrac{\tau}{2}\mathbb{E}\left[\textstyle\sum^{n}_{i=1}(x^{\prime}_{i,K}-\tfrac{1}{\tau}z_{i})^{2}\right]\geq\tfrac{\nu^{2}}{2\tau}\mathbb{E}\left[\mathbbold{1}_{x^{\prime}_{i,K}z_{i}<0}\right],$

where the expectation is taken over the quadratic functions of the form (4.9).
This means that we can construct a bound from below if we can get a grip on
the signs of each $z_{i}$. To that end, we follow the proof of [Sha13, Theorem
3]. The idea is to consider deterministic strategies that have only access to
a sequence of function evaluations. The KL-divergence will allow for relating
these function evaluations and the sign of $z_{i}$.

The key difference with respect to [Sha13], however, is the estimator. Given
some point $x_{k}$, our function evaluation $v_{k}$ is of the form
$v_{k}=\Im(f(x_{k}+i\delta u))+\xi$ for some $\delta>0$,
$u\in\mathbb{S}^{n-1}$ and noise realization $\xi$. Now observe that
$\Im(f_{z}(x+i\delta u))=\delta(\tau\langle x,u\rangle-\langle z,u\rangle)$.
Hence, conditioning on $z_{i}>0$ we get

$v_{k}=\delta\left(\tau\langle x_{k},u\rangle-\textstyle\sum_{j\neq
i}z_{j}u_{j}\right)-\nu u_{i}+\xi$ (4.10) whereas conditioning on $z_{i}<0$
yields $v_{k}=\delta\left(\tau\langle x_{k},u\rangle-\textstyle\sum_{j\neq
i}z_{j}u_{j}\right)+\nu u_{i}+\xi.$

Under the assumption that the noise is Gaussian one can now bound the KL-
divergence by $(2\nu u_{i})^{2}/(2\sigma_{\xi})$, e.g., see [Sha13, Lemma 5].
Using the fact that $u\in\mathbb{S}^{n-1}$ one can now exploit [Sha13, Lemma
4] and show that

$\tfrac{\nu^{2}}{2\tau}\mathbb{E}[\mathbbold{1}_{x^{\prime}_{i,K}z_{i}<0}]\geq\tfrac{n\nu^{2}}{4\tau}\left(1-\sqrt{\tfrac{2\nu^{2}K}{n\sigma_{\xi}}}\right).$

As such, selecting $\nu=\sqrt{({n\sigma_{\xi}})/({4K})}$ yields the desired
result. ∎

In the light of Theorem 4.6 and Remark 4.3, the following algorithms are rate
optimal. More specifically, one can show that the dependence on $\sigma_{\xi}$
is also optimal. Note that for quadratic functions we should not simply appeal
to Theorem 4.2 as that proof relies on $L_{2}(f)>0$.

######  Theorem 4.7 (Convergence rate of Algorithm 4.1 (b) with noise, $f$
being quadratic).

Let $f\in C^{\omega}(\mathcal{D})$ be a $\tau(f)$-strongly convex function
satisfying Assumption 3.1 for some $\bar{\delta}\in(0,1)$ and let
$\mathcal{K}\subset\mathcal{D}$ be a compact convex set. Suppose that $f$ has
a constant Hessian over $\mathcal{K}$, that is, (1.3) holds with $L_{2}(f)=0$.
Let $\\{x_{k}\\}_{k\geq 1}$ be the sequence of iterates generated by Algorithm
4.1 (b) with stepsize $\mu_{k}={2}/({\tau(f)k})$ and the sequence of smoothing
parameters defined for all $k\geq 1$ by $\delta_{k}=\delta$ with
$\delta\in(0,\kappa\bar{\delta}]$ for some $\kappa\in(0,1)$. Then, if the
oracle satisfies Assumption 1.2, the uniformly-averaged iterate
$\bar{x}_{K}=K^{-1}\sum^{K}_{k=1}x_{k}$ achieves the optimization error

$\mathbb{E}[f(\bar{x}_{K})-f(x^{\star})]\leq{O}\left(\frac{n}{\tau(f)}K^{-1}\right)+\widetilde{O}\left(\frac{n^{2}\sigma_{\xi}}{\tau(f)\delta^{2}}K^{-1}\right).$

###### Proof.

We can mainly follow the proof of Theorem 4.2, which relies itself largely on
[APT20]. To that end, let again $\sup_{x\in\mathcal{K}}\|\nabla f(x)\|_{2}\leq
G$ and set $a_{k}=\|x_{k}-x^{\star}\|_{2}^{2}$ such that

$\displaystyle\begin{array}[]{lcl}\mathbb{E}[f(x_{k})-f(x^{\star})|x_{k}]&\leq&\left\|\mathbb{E}[g_{\delta_{k}}(x_{k})|x_{k}]-\nabla
f(x_{k})\right\|_{2}\|x_{k}-x^{\star}\|_{2}+\tfrac{1}{2\mu_{k}}\mathbb{E}[a_{k}-a_{k+1}|x_{k}]\\\
&&\tfrac{\mu_{k}}{2}\mathbb{E}[\|g_{\delta_{k}}(x_{k})\|_{2}^{2}|x_{k}]-\tfrac{\tau(f)}{2}\mathbb{E}[a_{k}|x_{k}]\\\
&\overset{\eqref{equ:grad:delta:approx:error}}{\leq}&\tfrac{1}{2\mu_{k}}\mathbb{E}[a_{k}-a_{k+1}|x_{k}]+\tfrac{\mu_{k}}{2}\mathbb{E}[\|g_{\delta_{k}}(x_{k})\|_{2}^{2}|x_{k}]-\tfrac{\tau(f)}{2}\mathbb{E}[a_{k}|x_{k}].\end{array}$

Next, let $r_{k}=\mathbb{E}[a_{k}]$ such that we can write

$\mathbb{E}[f(x_{k})-f(x^{\star})]\leq\tfrac{1}{2\mu_{k}}(r_{k}-r_{k+1})-\tfrac{\tau(f)}{2}r_{k}+\tfrac{\mu_{k}}{2}\mathbb{E}[\|g_{\delta_{k}}(x_{k})\|_{2}^{2}].$
(4.11)

To allow for an identical stepsize as before, we replace $-\tau(f)/2$ with
$-\tau(f)/4$. Summing (4.11) over $k$ yields

$\displaystyle\textstyle\sum^{K}_{k=1}\mathbb{E}[f(x_{k})-f(x^{\star})]\leq$
$\displaystyle\tfrac{1}{2}\textstyle\sum^{K}_{k=1}\left(\tfrac{1}{\mu_{k}}(r_{k}-r_{k+1})-\tfrac{\tau(f)}{2}r_{k}\right)+\textstyle\sum^{K}_{k=1}\tfrac{\mu_{k}}{2}\mathbb{E}[\|g_{\delta_{k}}(x_{k})\|_{2}^{2}].$

As we selected $\mu_{k}=2/(\tau(f)k)$ we can simplify the above by using the
same reasoning as in [APT20], that is

$\displaystyle\textstyle\sum^{K}_{k=1}\left(\tfrac{1}{\mu_{k}}(r_{k}-r_{k+1})-\tfrac{\tau(f)}{2}r_{k}\right)\leq
r_{1}\left(\tfrac{1}{\mu_{1}}-\tfrac{\tau(f)}{2}\right)+\textstyle\sum^{K}_{k=2}r_{k}\left(\tfrac{1}{\mu_{k}}-\tfrac{1}{\mu_{k-1}}-\tfrac{\tau(f)}{2}\right)=0.$

Indeed, without the scaling of $\tau(f)$ our stepsize would have been
$\mu_{k}=1/(\tau(f)k)$. Note that we rely on the $\tau(f)$-strong convexity.
Using the observation from above and plugging in the stepsize $\mu_{k}$
elsewhere yields

$\displaystyle\begin{array}[]{lcl}\textstyle\sum^{K}_{k=1}\mathbb{E}[f(x_{k})-f(x^{\star})]&\leq&\tfrac{1}{\tau(f)}\textstyle\sum^{K}_{k=1}\tfrac{1}{k}\mathbb{E}[\|g_{\delta_{k}}(x_{k})\|_{2}^{2}]\\\
&\leq&\tfrac{n^{2}}{\tau(f)}\textstyle\sum^{K}_{k=1}\tfrac{1}{k}\left[\tfrac{1}{n}\|\nabla
f(x_{k})\|_{2}^{2}+\frac{1}{\delta_{k}^{2}}\sigma_{\xi}\right].\end{array}$

Now, minimizing over $\\{\delta_{k}\\}_{k}$ clearly yields a desire to pick a
larger and fixed $\delta$ cf. Theorem 4.2. Combining this with the bound on
$\nabla f(x)$ yields by (A.1)

$\displaystyle\begin{array}[]{lcl}\textstyle\sum^{K}_{k=1}\mathbb{E}[f(x_{k})-f(x^{\star})]&\leq&\tfrac{n^{2}}{\tau(f)}\left[\tfrac{1}{n}G^{2}+\tfrac{1}{\delta^{2}}\sigma_{\xi}\right](\log(K)+1),\end{array}$

as such we obtain the optimization error
$\mathbb{E}[f(\bar{x}_{k})-f(x^{\star})]\leq\tfrac{n^{2}}{\tau(f)K}\left[\tfrac{1}{n}G^{2}+\tfrac{1}{\delta^{2}}\sigma_{\xi}\right](\log(K)+1)$.
∎

######  Example 4.8 (Numerical strongly-convex optimization).

Here we exemplify what can go wrong and how the proposed complex-step method
handles this. Consider for $n=10^{3}$ the problem of solving

$\displaystyle\operatorname*{\text{minimize}}_{x\in\mathbb{B}^{n}}$
$\displaystyle\tfrac{1}{2}\langle x,x\rangle.$

We let $\sigma_{\xi}=\epsilon_{M}^{4}$, with $\delta=1$ or $\delta=10^{-100}$
(two extremes) and compare Theorem 4.7 (CS algorithm) against a state-of-the-
art multi-point method [APT20, Theorem 5.1] ($\beta$ algorithm)121212The small
$\sigma_{\xi}$ aides the exposition as a larger $\sigma_{\xi}$ would merely
delay the effect.. Their stepsize equals ours, but their smoothing parameter
equals

$\delta_{k}=\left(\frac{3n^{2}\sigma_{\xi}}{4k+9n^{2}}\right)^{\tfrac{1}{4}},\quad
k=1,2,\dots,K.$

In Figure 1(b) we show for 250 experiments
($x_{1}\overset{i.i.d.}{\sim}\mathcal{N}(0,I_{n})$), the differences in
convergence. Indeed, the proposed complex-step does not suffer from
cancellation errors as can be seen in Figure D.1. The reason why
$\delta=10^{-100}$ works (unreasonably) well is due to the constraints and the
averaging, each iteration lives on $\partial\mathbb{B}^{n}$ (recall error
terms of the form $1/\delta^{2}$). Although the setting is somewhat esoteric,
this example does show the possibility of catastrophic cancellation and how to
resolve it.

Similar to Theorem 4.4, we analyze unconstrained zeroth-order optimization
when $f$ is quadratic.

######  Theorem 4.9 (Convergence rate of Algorithm 4.1 (a) with noise, $f$
being quadratic).

Let $f\in C^{\omega}(\mathcal{D})$ be a $\tau(f)$-strongly convex quadratic
function satisfying Assumption 3.1 for some $\bar{\delta}\in(0,1)$ with
$x^{\star}\in\mathrm{int}(\mathcal{D})$. Suppose that $f$ has a Lipschitz
gradient and constant Hessian, that is, (1.2) and (1.3) hold, for $L_{1}(f)>0$
and $L_{2}(f)=0$, respectively. Let $\\{x_{k}\\}_{k\geq 1}$ be the sequence of
iterates generated by Algorithm 4.1 (a) for

$\displaystyle\begin{array}[]{llll}&\mu_{k}=\tfrac{1}{\tau(f)K},&\delta_{k}=\delta,&k=1,\dots,K_{0},\\\
&\mu_{k}=\tfrac{2}{\tau(f)k},&\delta_{k}=\delta,&k=K_{0}+1,\dots,K,\end{array}$

with $K_{0}=\left\lfloor\tfrac{4nL_{1}(f)^{2}}{\tau(f)^{2}}\right\rfloor$ and
$\delta\in(0,\kappa\bar{\delta}]$ for some $\kappa\in(0,1)$. Then, if the
oracle satisfies Assumption 1.2 and $K\geq 2K_{0}$ we incur for
$\bar{x}_{K_{0},K}=\tfrac{1}{K-K_{0}}\sum^{K}_{k=K_{0}+1}x_{k}$ the
optimization error

$\displaystyle\mathbb{E}[f(\bar{x}_{K_{0},K})-f(x^{\star})]$
$\displaystyle\leq
O\left(\frac{nL_{1}(f)^{2}}{\tau(f)}\|x_{1}-x^{\star}\|_{2}^{2}K^{-1}\right)+\widetilde{O}\left(\frac{n^{2}\sigma_{\xi}}{\tau(f)\delta^{2}}K^{-1}\right).$
(4.12)

###### Proof.

The proof will be a combination of Theorem 4.4 and Theorem 4.7.

Again, set $a_{k}=\|x_{k}-x^{\star}\|_{2}^{2}$, then

$\displaystyle\begin{array}[]{lcl}\mathbb{E}[f(x_{k})-f(x^{\star})|x_{k}]&\overset{\eqref{equ:grad:delta:approx:error}}{\leq}&\tfrac{1}{2\mu_{k}}\mathbb{E}[a_{k}-a_{k+1}|x_{k}]+\tfrac{\mu_{k}}{2}\mathbb{E}[\|g_{\delta_{k}}(x_{k})\|_{2}^{2}|x_{k}]-\tfrac{\tau(f)}{2}\mathbb{E}[a_{k}|x_{k}].\end{array}$

Next, let $r_{k}=\mathbb{E}[a_{k}]$ such that by $\|\nabla
f(x_{k})\|_{2}^{2}\leq L_{1}(f)^{2}\|x_{k}-x^{\star}\|_{2}^{2}$ we can write

$\displaystyle\mathbb{E}[f(x_{k})-f(x^{\star})]$
$\displaystyle\leq\tfrac{1}{\mu_{k}}(r_{k}-r_{k+1})-{\tau(f)}r_{k}+{\mu_{k}}\left(nL_{1}(f)^{2}r_{k}+\tfrac{n^{2}}{\delta_{k}^{2}}\sigma_{\xi}\right)$
(4.13)

Now we use the step- and smoothingsize for $k=1,\dots,K_{0}$, that is,
$\mu_{k}=1/(\tau(f)K)$, $\delta_{k}=\delta$, and observe that

$\displaystyle r_{k+1}$ $\displaystyle\leq
r_{k}-\tau(f)\mu_{k}r_{k}+\mu_{k}^{2}\left(nL_{1}(f)^{2}r_{k}+\tfrac{n^{2}}{\delta_{k}^{2}}\sigma_{\xi}\right)=\left(1-\tfrac{1}{K}+\tfrac{L_{1}(f)^{2}}{(\tau(f)K)^{2}}n\right)r_{k}+\nu_{K}=a_{K}r_{k}+\nu_{K}$

for $a_{k}$ as between brackets and $\nu_{K}$ defined as

$\displaystyle\nu_{K}=$
$\displaystyle\tfrac{n^{2}}{(\delta\tau(f)K)^{2}}\sigma_{\xi}.$

We now proceed with bounding $r_{K_{0}+1}$. As in [APT20], set

$q_{K}=1+\tfrac{L_{1}(f)^{2}}{(\tau(f)K)^{2}}n$

by iterating over $r_{k}$ it follows that

$\displaystyle r_{K_{0}+1}\leq
a_{K}^{K_{0}}r_{1}+\textstyle\sum^{K_{0}-1}_{i=0}a_{K}^{i}\nu_{K}\leq\left(r_{1}+\tfrac{(\tau(f)K)^{2}}{L_{1}(f)^{2}n}\nu_{K}\right)q_{K}^{K_{0}}.$

Now assume that that $K_{0}$ is as in the theorem, then as $\log(1+x)\leq x$

$\displaystyle q_{K}^{K_{0}}=$
$\displaystyle\mathrm{exp}\left(K_{0}\log\left(1+\tfrac{L_{1}(f)^{2}}{(\tau(f)K)^{2}}n\right)\right)\leq\mathrm{exp}\left(\tfrac{4nL_{1}(f)^{2}}{\tau(f)^{2}}\log\left(1+\tfrac{L_{1}(f)^{2}}{(\tau(f)K)^{2}}n\right)\right)\leq\mathrm{exp}\left(\tfrac{4n^{2}L_{1}(f)^{4}}{\tau(f)^{4}K^{2}}\right).$

Fix any $\bar{C}\in(0,\tfrac{1}{16})$, when

$K\geq\sqrt{\tfrac{4n^{2}L_{1}(f)^{4}}{\tau(f)^{4}\bar{C}}}$

then $K\geq 2K_{0}$ and $q^{K_{0}}_{K}\leq e^{\bar{C}}=C_{4}$. As such,

$\displaystyle r_{K_{0}+1}\leq$
$\displaystyle\left(r_{1}+\tfrac{(\tau(f)K)^{2}}{L_{1}(f)^{2}n}\nu_{K}\right)C_{4}\leq\left(r_{1}+\tfrac{(\tau(f)K)^{2}}{L_{1}(f)^{2}n}\tfrac{n^{2}}{(\delta\tau(f)K)^{2}}\sigma_{\xi}\right)C_{4}=\left(r_{1}+\tfrac{n}{\delta^{2}L_{1}(f)^{2}}\sigma_{\xi}\right)C_{4}.$

Now we return to our normal step- and smoothingsizes, that is
$\mu_{k}=2/(\tau(f)k)$, $\delta_{k}=\delta$, for $k\geq K_{0}+1$. By plugging
this into (4.13) we get

$\displaystyle(K-K_{0})\mathbb{E}[f(\bar{x}_{K_{0},K})-f(x^{\star})]\leq$
$\displaystyle\textstyle\sum^{K}_{k=K_{0}+1}\tfrac{\tau(f)k}{2}(r_{k}-r_{k+1})-{\tau(f)}r_{k}+\tfrac{2}{\tau(f)k}nL_{1}(f)^{2}r_{k}$
$\displaystyle+$
$\displaystyle\textstyle\sum^{K}_{k=K_{0}+1}\tfrac{2}{\tau(f)k}\tfrac{n^{2}}{\delta^{2}}\sigma_{\xi}$

By construction of $K_{0}$ we have that for $k\geq K_{0}+1$,
$\tau(f)/2\geq(2nL_{1}(f)^{2})/(\tau(f)k)$. Hence

$\displaystyle(K-K_{0})\mathbb{E}[f(\bar{x}_{K_{0},K})-f(x^{\star})]$
$\displaystyle\leq\tfrac{\tau(f)}{2}\textstyle\sum^{K}_{k=K_{0}+1}k(r_{k}-r_{k+1})-r_{k}+U_{K_{0},K}.$

where by (A.1)

$\displaystyle
U_{K_{0},K}=\textstyle\sum^{K}_{k=K_{0}+1}\tfrac{2}{\tau(f)k}\tfrac{n^{2}}{\delta^{2}}\sigma_{\xi}\leq\tfrac{2}{\tau(f)}\tfrac{n^{2}}{\delta^{2}}\sigma_{\xi}(\log(K)+1).$

As demonstrated in [APT20], one can now construct the bound
$\textstyle\sum^{K}_{k=K_{0}+1}k(r_{k}-r_{k+1})-r_{k}\leq K_{0}r_{K_{0}+1}$
where the last term is exactly the term we could bound before. In combination
with the bound on $K_{0}$ itself, we find that

$\displaystyle(K-K_{0})\mathbb{E}[f(\bar{x}_{K_{0},K})-f(x^{\star})]$
$\displaystyle\leq\tfrac{\tau(f)}{2}\tfrac{4nL_{1}(f)^{2}}{\tau(f)^{2}}\left(r_{1}+\tfrac{n}{\delta^{2}L_{1}(f)^{2}}\sigma_{\xi}\right)C_{4}+\tfrac{2}{\tau(f)}\tfrac{n^{2}}{\delta^{2}}\sigma_{\xi}(\log(K)+1).$

By our selection of $\bar{C}$ we have that $K\geq 2K_{0}$ and as such

$\displaystyle\mathbb{E}[f(\bar{x}_{K_{0},K})-f(x^{\star})]$
$\displaystyle\leq\tfrac{2nL_{1}(f)^{2}}{\tau(f)K}\left(r_{1}+\tfrac{n}{\delta^{2}L_{1}(f)^{2}}\sigma_{\xi}\right)C_{4}+\tfrac{2n^{2}}{\delta^{2}\tau(f)K}\sigma_{\xi}(\log(K)+1).$

Now, reordering terms yields (4.12). ∎

It is important to highlight that the stepsizes for the quadratic cases are
identical to the general case. As such, no knowledge of the quadratic nature
is required, but if $f$ happens to be quadratic, the algorithm performs
optimally.

###  4.4 Online optimization

Online optimization shows up in settings where the objective might change due
to the presence of more information, say, when more data becomes available. In
the online case one is interested in bounding the regret of the form

$\tfrac{1}{K}\textstyle\sum^{K}_{k=1}\mathbb{E}[f_{k}(x_{k})]-\inf_{x\in\mathcal{D}}\tfrac{1}{K}\sum^{K}_{k=1}f_{k}(x).$
(4.14)

As in for example [BP16], the proof techniques are largely the same as for the
stochastic cases above. We consider the following setting to exemplify the
possibilities. Note, here the algorithm proceeds as

$x_{k+1}=\Pi_{\mathcal{K}}\left(x_{k}-\mu_{k}\cdot\left(\tfrac{n}{\delta_{k}}\Im\left(f_{k+1}(x_{k}+i\delta_{k}u_{k})\right)u_{k}+\tfrac{n}{\delta_{k}}\xi_{k}u_{k}\right)\right).$

######  Theorem 4.10 (Online optimization, convergence rate of Algorithm 4.1
(b) with noise, $f_{k}$ being quadratic).

Let all $f_{k}\in C^{\omega}(\mathcal{D})$ be $\tau(f_{k})$-strongly convex
functions satisfying Assumption 3.1 for some $\bar{\delta}\in(0,1)$ and let
$\mathcal{K}\subset\mathcal{D}$ be a compact convex set. Suppose that all
$f_{k}$ have a mutual Lipschitz gradient and a constant Hessian over
$\mathcal{K}$, that is, (1.2) and (1.3) hold, for some constant
$L_{1}(f_{k})>0$ and $L_{2}(f_{k})=0$, respectively. Set
$\tau=\min_{k}\tau(f_{k})$ and let $\\{x_{k}\\}_{k\geq 1}$ be the sequence of
iterates generated by Algorithm 4.1 (b) with stepsize $\mu_{k}={2}/({\tau k})$
and the sequence of smoothing parameters defined for all $k\geq 1$ by
$\delta_{k}=\delta$ with $\delta\in(0,\kappa\bar{\delta}]$ for some
$\kappa\in(0,1)$. Then, if the oracle satisfies Assumption 1.2 we incur the
regret

$\tfrac{1}{K}\textstyle\sum^{K}_{k=1}\mathbb{E}[f_{k}(x_{k})]-\inf_{x\in\mathcal{K}}\tfrac{1}{K}\sum^{K}_{k=1}f_{k}(x)\leq\widetilde{O}\left(\frac{n^{2}}{\tau}K^{-1}\right)$
(4.15)

###### Proof.

As (4.5) holds for any $x$ and not just $x^{\star}$, the proof is effectively
identical to that of Theorem 4.7. By appealing to that proof we have the
following bound immediately

$\displaystyle\begin{array}[]{lcl}\textstyle\sum^{K}_{k=1}\mathbb{E}[f_{k}(x_{k})-f_{k}(x)]&\leq&\tfrac{n^{2}}{\tau}\left[\tfrac{1}{n}G^{2}+\tfrac{1}{\delta^{2}}\sigma_{\xi}\right](\log(K)+1)\end{array}$

and as such we obtain the regret bound (4.15). ∎

###  4.5 Numerical estimation of $\tau(f)$

As most regret bounds and stepsizes contain terms of the form $1/\tau(f)$ one
should take care in estimating the strong convexity parameter $\tau(f)$. An
arbitrarily small $\tau(f)$ complies with the definition but could lead for
instance to numerical overflow due to large stepsizes. In fact, it is known
that either under- or overestimating $\tau(f)$ can have detrimental effects on
convergence properties, especially in accelerated schemes [OC15]. When one has
access to gradients, line-search-like schemes are possible to estimate both
$\tau(f)$ and $L_{1}(f)$ [Nes13]. However, when the gradient direction is
random, this is less straight-forward.

Fitting a quadratic model using a (recursive) least squares approach can
grossly overestimate $\tau(f)$. For example, consider the function
$f(x)=x^{4}+x^{2}+\tfrac{1}{2}\lambda x^{2}$. One might have access to
$\lambda>0$, e.g., by means of being a regularization parameter. Then, fitting
a quadratic model to this function yields (asymptotically) a strong convexity
estimate of $6+\lambda$ instead of $2+\lambda$.

We propose simple routines to estimate the largest $\tau(f)$ satisfying (4.1),
denoted $\bar{\tau}(f)$. Here we exploit the fact that we have a sequence of
function evaluations, which remain commonly and unfortunately unused in this
line of zeroth-order optimization schemes. We also assume to have knowledge of
some lower bound $\tau_{0}>0$ such that $\tau_{0}\leq\bar{\tau}(f)$, which is
frequently available due to regularization. In terms of the dimension $n$ we
identify two regimes, small (medium) scale $n\leq 10^{3}$ and large scale
$n>10^{3}$.

  1. (i)

(Small scale): Using the data at hand we can construct an explicit quadratic
model in $(P,q,r)$ that bounds $f(x)$ from below. Due to the inherent
randomness, $f\in C^{\omega}$, and the possibility of selecting $\delta_{k}$
close to $0$, one has (for $\delta_{k}$ sufficiently small) a sufficiently
accurate quadratic model by using $N(n)\geq\tfrac{1}{2}(n+1)(n+2)$ data-points
in the following semidefinite program (SDP)

$\displaystyle\operatorname*{\text{minimize}}_{(P,q,r)}$
$\displaystyle\textstyle\sum^{K^{\prime}+N(n)-1}_{k=K^{\prime}}\Re(f(x_{k}+i\delta_{k}u_{k}))-\langle\tfrac{1}{2}Px_{k}+q,x_{k}\rangle-r+\delta_{k}^{2}\langle\tfrac{1}{2}Pu_{k},u_{k}\rangle$
(4.16) $\displaystyle\operatorname*{\text{subject to}}$ $\displaystyle
P\in\mathcal{S}^{n}_{\succ
0},\,q\in\mathbb{R}^{n},\,r\in\mathbb{R},\,P\succeq\tau_{0}I_{n},$
$\displaystyle\Re(f(x_{k}+i\delta_{k}u_{k}))\geq\langle\tfrac{1}{2}Px_{k}+q,x_{k}\rangle+r-\delta_{k}^{2}\langle\tfrac{1}{2}Pu_{k},u_{k}\rangle,$
$\displaystyle\text{for }k=K^{\prime},\dots,K^{\prime}+N(n)-1,$

for some $K^{\prime}\geq 1$. Now, an approximation of $\bar{\tau}(f)$ follows
by setting $\widehat{\tau}(f)=\min_{i}\\{\lambda_{i}(P^{\star})\\}$. Indeed,
models as such can now also be used to further fine-tune the proposed
algorithms.

Let us elaborate on the aforementioned claims in the unconstrained case. The
constrained case is less predictable. First of all, we want to have a tight
quadratic model that approximates the function $f$ from below. To use our data
economically, we have to do with samples of the form
$\Re(f(x_{k}+i\delta_{k}u_{k}))$ instead of $f(x_{k})$. Note, these samples
might be corrupted by noise. Then, the objective in combination with the
inequality constraints in (4.16) enforce that $(P,q,r)$ parametrizes a
quadratic model, approximately from below, that is as close as possible to the
available data. To parametrize this model one needs at most
$n(n+1)/2+n+1=\tfrac{1}{2}(n+1)(n+2)$ data-points indeed. Here, a data-point
compromises the $4$-tuple
$(x_{k},u_{k},\delta_{k},f(x_{k}+i\delta_{k}u_{k}))$. Now, as we sample
$u_{1},u_{2},\dots$ uniformly and independently from $\mathbb{S}^{n-1}$, the
set $\\{u_{1},u_{2},\dots,u_{n}\\}$ will $\mathbb{P}$-a.s. span
$\mathbb{R}^{n}$. Then, as the noise terms $\xi_{1},\xi_{2},\dots$ that
potentially enter the oracle are independent of $u_{1},u_{2},\dots$ and $f\neq
0$ can only vanish on sets of measure $0$ by the real-analytic assumption, we
must have that $x_{k}$ is $\mathbb{P}$-a.s. not parallel to $x_{k+1}$, cf.
Algorithm 4.1.

  2. (ii)

(Large scale): When $n$ is large, we follow the ideas as set forth in [AM19].
Denote by $\mathsf{dd}^{n}$ the diagonally dominant matrices in
$\mathbb{R}^{n\times n}$. That is, $A\in\mathsf{dd}^{n}$ when
$a_{ii}\geq\sum_{j\neq i}|a_{ij}|$ for all $i\in[n]$. This allows for a
polytopic representation of the constraint $P-\tau_{0}I_{n}\succeq 0$. Now we
transform (4.16) in the diagonally dominant program (DDP) by identifying
$p\in\mathbb{R}^{n(n+1)/2}$ with $\mathrm{svec}(P)$, that is, $P$ is not an
additional decision variable but merely an auxiliary variable to simplify
notation

$\displaystyle\operatorname*{\text{minimize}}_{(p,q,r)}$
$\displaystyle\textstyle\sum^{K^{\prime}+N(n)-1}_{k=K^{\prime}}\Re(f(x_{k}+i\delta_{k}u_{k}))-\tfrac{1}{2}\langle
p,x_{k}\otimes_{s}x_{k}-\delta_{k}^{2}u_{k}\otimes_{s}u_{k}\rangle-\langle
q,x_{k}\rangle-r$ (4.17) $\displaystyle\operatorname*{\text{subject to}}$
$\displaystyle
p\in\mathbb{R}^{n(n+1)/2},\,q\in\mathbb{R}^{n},\,r\in\mathbb{R},\,P=\mathrm{smat}(p),$
$\displaystyle P_{ii}-\tau_{0}\geq\textstyle\sum_{j\neq i}|P_{ij}|,\,i=1,\dots
n,$ $\displaystyle\Re(f(x_{k}+i\delta_{k}u_{k}))\geq\tfrac{1}{2}\langle
p,x_{k}\otimes_{s}x_{k}-\delta_{k}^{2}u_{k}\otimes_{s}u_{k}\rangle+\langle
q,x_{k}\rangle+r,$ $\displaystyle\text{for
}k=K^{\prime},\dots,K^{\prime}+N(n)-1.$

Here, $\otimes_{s}$ denotes the symmetric Kronecker product. See [MHA20] for a
recent survey on large-scale SDPs and Section B for more on the DDP-based
approximation of $\bar{\tau}(f)$. Specifically, we can iteratively improve the
basis in (4.17), such that $p^{\star}$ with respect to (4.17) converges weakly
to $P^{\star}$ with respect to (4.16).

######  Example 4.11 (Numerical performance of $\tau(f)$ estimation).

To show how the proposed estimation scheme for $\tau(f)$ can be beneficial we
look at a transparent (closed-form solutions are available) example. Consider
the $\ell_{2}$-regularized least-squares problem

$\operatorname*{\text{minimize}}_{x\in
r\mathbb{B}^{n}}\quad\tfrac{1}{2}\langle
Ax-b,Ax-b\rangle+\tfrac{1}{2}\lambda\langle x,x\rangle$

for $r>0$ such that $x^{\star}\in\mathrm{int}(r\mathbb{B}^{n})$. In many
problems one might have knowledge of the regularization parameter $\lambda>0$
but not of the remaining objective terms. As such, we start with
$\tau_{0}=\lambda$ and use the SDP formulation (4.16) to approximate
$\bar{\tau}(f)$ from below by $\widehat{\tau}(f)$. We do 250 experiments
($x_{1}\overset{i.i.d.}{\sim}\mathcal{N}(0,I_{n})$,
$\mathrm{vec}(A)\overset{i.i.d.}{\sim}\mathcal{N}(0,I_{mn})$,
$b\overset{i.i.d.}{\sim}\mathcal{N}(0,I_{m})$) for $n=10$, $m=20$ and
$\lambda=10^{-4}$. We plug the estimation scheme into Algorithm 4.1 (b) (for
$\delta_{k}=\epsilon_{M}$ and $\sigma_{\xi}=\epsilon_{M}^{4}$), that is,
compute $\widehat{\tau}(f)$ once at $K=\tfrac{1}{2}(n+1)(n+2)$, and show the
results in Figure 1(c). The approximation clearly speeds up the convergence
and closely resembles that under $\bar{\tau}(f)$. Section B (Appendix)
presents a similar example for (4.17).

The take away of this section is not only a routine to estimate $\tau(f)$, but
also the observation that this can be done directly using the complex function
evaluations of the form $f(x_{k}+i\delta_{k}u_{k})$.

##  5 Outlook: nonconvex zeroth-order optimization

At last we consider a critical point in a possibly non-convex program. Note,
we do not assume that our function is locally convex. We exploit that the
gradient of $f\in C^{\omega}(\mathcal{D})$ is uniformly bounded over any
compact subset of $\mathcal{D}$.

######  Theorem 5.1 (Convergence rate of Algorithm 4.1 (b) to a critical
point).

Let $f\in C^{\omega}(\mathcal{D})$ be a — not necessarily convex — function
that satisfies Assumption 3.1 for some $\bar{\delta}\in(0,1)$. Suppose that
$f$ has a Lipschitz gradient and Hessian on $\mathcal{K}\subset\mathcal{D}$,
that is, (1.2) and (1.3) hold, for some constants $L_{1}(f)>0$ and
$L_{2}(f)\geq 0$, respectively. Let $\\{x_{k}\\}_{k\geq 1}$ be the sequence of
iterates generated by Algorithm 4.1 (b) with stepsize
$\mu_{k}={1}/({nL_{1}(f)k^{2/3}})$ and the sequence of smoothing parameters
defined for all $k\geq 1$ by $\delta_{k}=\delta k^{-1/6}$ with
$\delta\in(0,\kappa\bar{\delta}]$ for some $\kappa\in(0,1)$. Let $x^{\star}$
be a global minimum of $f$, then,

$\displaystyle\min_{k\geq 1}\mathbb{E}[\|\nabla f(x_{k})\|^{2}_{2}]\leq
O\left({nL_{1}(f)(f(x_{1})-f(x^{\star}))}K^{-\tfrac{1}{3}}\right)+\widetilde{O}\left(n{(\delta^{2}+\frac{\sigma_{\xi}}{\delta^{2}})}K^{-\tfrac{1}{3}}\right).$
(5.1)

###### Proof.

Our proof will be similar to constructions as set forth in [Nes03]. As $f\in
C^{1,1}_{L_{1}(f)}$ one has

$\displaystyle\begin{array}[]{lcl}f(x_{k+1})&\overset{\eqref{equ:L1:upper}}{\leq}&f(x_{k})-\mu_{k}\langle\nabla
f(x_{k}),g_{\delta_{k}}(x_{k})\rangle+\tfrac{1}{2}\mu_{k}^{2}L_{1}(f)\|g_{\delta_{k}}(x_{k})\|_{2}^{2}\\\
&=&f(x_{k})-\mu_{k}\|\nabla f(x_{k})\|_{2}^{2}-\mu_{k}\langle\nabla
f(x_{k}),g_{\delta_{k}}(x_{k})-\nabla
f(x_{k})\rangle+\tfrac{1}{2}\mu_{k}^{2}L_{1}(f)\|g_{\delta_{k}}(x_{k})\|_{2}^{2}.\end{array}$

Now taking expectation, applying the Cauchy-Schwarz inequality and using both
(3.7) and (3.8) results in

$\displaystyle\mathbb{E}_{u_{k}\sim\mathbb{S}^{n-1}}\left[f(x_{k+1})|x_{k}\right]\leq$
$\displaystyle f(x_{k})-\mu_{k}\|\nabla
f(x_{k})\|_{2}^{2}+\mu_{k}{C_{1}n\delta_{k}^{2}}\|\nabla f(x_{k})\|_{2}$
$\displaystyle+\tfrac{1}{2}\mu_{k}^{2}L_{1}(f)\left(n\|\nabla
f(x_{k})\|_{2}^{2}+{C_{2}n^{2}\delta_{k}^{4}}+C_{3}n^{2}\delta_{k}^{2}\|\nabla
f(x_{k})\|_{2}+\tfrac{n^{2}}{\delta_{k}^{2}}\sigma_{\xi}\right).$

Then, taking expectation over $u_{1},\dots,u_{k-1}$, plugging in our stepsize
$\mu_{k}=1/(nL_{1}(f)k^{\tfrac{2}{3}})$ applying Jensen’s inequality and
rearranging yields

$\displaystyle k^{-\tfrac{2}{3}}\tfrac{1}{2nL_{1}(f)}\mathbb{E}\left[\|\nabla
f(x_{k})\|_{2}^{2}\right]\leq$ $\displaystyle\mathbb{E}[f(x_{k})-f(x_{k+1})]$
(5.2)
$\displaystyle+k^{-\tfrac{2}{3}}\tfrac{C_{4}\delta_{k}^{2}}{L_{1}(f)}\mathbb{E}\left[\|\nabla
f(x_{k})\|_{2}\right]+k^{-\tfrac{4}{3}}\tfrac{C_{2}\delta_{k}^{4}}{2L_{1}(f)}+k^{-\tfrac{4}{3}}\tfrac{1}{2L_{1}(f)\delta_{k}^{2}}\sigma_{\xi}.$

As we consider a global minimum we have that $f(x_{k})\geq f(x^{\star})$. Now,
define $\phi_{k}=(1/k^{\tfrac{2}{3}})\mathbb{E}[\|\nabla f(x_{k})\|^{2}_{2}]$
and $\sup_{x\in\mathcal{K}}\|\nabla f(x)\|_{2}=G<+\infty$, then, a telescoping
argument yields

$\displaystyle\textstyle\sum^{K}_{k=1}\phi_{K}\leq
2nL_{1}(f)(f(x_{1})-f(x^{\star}))+2nC_{4}G\sum^{K}_{k=1}k^{-\tfrac{2}{3}}\delta_{k}^{2}+nC_{2}\sum^{K}_{k=1}k^{-\tfrac{4}{3}}\delta_{k}^{4}+\sum^{K}_{k=1}k^{-\tfrac{4}{3}}\tfrac{n}{\delta_{k}^{2}}\sigma_{\xi}.$

Now plug in $\delta_{k}=\delta k^{-\tfrac{1}{6}}$ and get

$\displaystyle\textstyle\sum^{K}_{k=1}\phi_{K}$ $\displaystyle\leq\textstyle
2nL_{1}(f)(f(x_{1})-f(x^{\star}))+{2nC_{4}G}\sum^{K}_{k=1}k^{-1}\delta^{2}+{nC_{2}}\sum^{K}_{k=1}k^{-2}\delta^{4}+\sum^{K}_{k=1}k^{-1}\tfrac{n}{\delta^{2}}\sigma_{\xi}$
$\displaystyle\leq
2nL_{1}(f)(f(x_{1})-f(x^{\star}))+C_{5}(\log(K)+1)(nG\delta^{2}+\tfrac{n}{\delta^{2}}\sigma_{\xi})+C_{6}n\delta^{4}.$

As such $\sum^{K}_{k=1}\phi_{K}\leq h(K)$, for $h(K)$ corresponding to the
right-most term above. This implies that $\min_{k}\phi_{k}\leq{h(K)/K}$. By
definition of $\phi_{k}$ we have $K^{-\tfrac{2}{3}}\min_{k}\mathbb{E}[\|\nabla
f(x_{k})\|^{2}_{2}]\leq\min_{k}\phi_{k}$. Combing these observations yields

$\min_{k\geq 1}\mathbb{E}[\|\nabla f(x_{k})\|^{2}_{2}]\leq
K^{-\tfrac{1}{3}}\left(2nL_{1}(f)(f(x_{1})-f(x^{\star}))+C_{5}n(\log(K)+1)(G\delta^{2}+\tfrac{\sigma_{\xi}}{\delta^{2}})+C_{6}n\delta^{4}\right).$

∎

Although the rate is relatively slow cf. [GL13], the approach appears to be
scalable, in contrast to common Monte Carlo methods [PS17]. Sharpening and
further generalizing Theorem 5.1 is left for future work.

Now we provide an numerical experiment, showing that Theorem 5.1 can handle
the noise, in contrast to the nonconvex algorithm as proposed in [JYK21] for
the deterministic setting.

######  Example 5.2 (Himmelblau function).

Consider optimizing a Himmelblau function over a closed ball centred at $0$,
in particular, consider

$\operatorname*{\text{minimize}}_{x\in
6\mathbb{B}^{2}}\quad\left(\big{(}x^{(1)}\big{)}^{2}+x^{(2)}-11\right)^{2}+\left(x^{(1)}+\big{(}x^{(2)}\big{)}^{2}-7\right)^{2}.$
(5.3)

The minimum value of (5.3) is $f^{\star}=0$. We will compare [JYK21, Algorithm
1] with stepsize $\mu_{k}=\mu=1/(nL_{1}(f))$ against Theorem 5.1. That is, we
compare a plain nonconvex algorithm (with $\delta_{k}=10^{-6}/k)$ against its
counterpart that is designed to handle noise (with
$\delta_{k}=10^{-6}k^{-1/6}$). We consider $8$ initial conditions (circles)
and show the results in Figure 5.1. The dark stars indicate minima of $f$,
whereas the light star is merely a local minima. We see that the algorithm
adapted to the noise can handle the perturbations well whereas the other
algorithm diverges. Note that formalizing these observations is left for
future work.

(a) The error $f(\bar{x}_{K})-f^{\star}$

(b) The trajectories $\bar{x}_{1},\dots,\bar{x}_{K}$.

Figure 5.1: Numerical outcomes of Example 5.2.

We end with an example pertaining to partial differential equations (PDEs).
PDEs are relevant as on the one hand, closed-form solutions are rare and
numerical solutions (approximations) are often a necessity, on the other hand,
analyticity of solutions has been studied since the early 50s, see for example
[Mor58, Mor58a].

######  Example 5.3 (PDE-constrained optimization).

PDEs can rarely be solved in closed-form and one commonly resorts to numerical
schemes, however, schemes that often lend themselves to the complex-lifting as
set forth in this article. In this example we show that there are already
examples that meet the conditions of Theorem 5.1. In particular, let $u$ be a
velocity field on $\mathbb{R}^{2}$, with abuse of notation, $(x,y)$ denote the
usual coordinates on $\mathbb{R}^{2}$. This velocity field is induced by a
solid sphere in $\mathbb{R}^{2}$ that moves in the negative $x$-direction with
a velocity $V$. We are interested in finding the optimal radius $r$ of this
sphere such that norm of the velocity field at the point
$(2,2)\in\mathbb{R}^{2}$ is minimized. When constraining the radius to the
interval $[1,2]$, then under idealized conditions (incompressibility and
irrotationality), we can consider the following PDE-constrained optimization
problem

$\displaystyle\operatorname*{\text{minimize}}_{r,u}$
$\displaystyle\|u(2,2)\|_{2}^{2}$ (5.4)
$\displaystyle\operatorname*{\text{subject to}}$
$\displaystyle\mathrm{div}(u)=0,\,\mathrm{curl}(u)=0,$ $\displaystyle\langle
u,n_{x}\rangle=\langle V,n_{x}\rangle,\,\forall x\in
r\mathbb{S}^{1},\,n_{x}\in(T_{x}r\mathbb{S}^{1})^{\perp},$ $\displaystyle
r\in[1,2],\,u\in C^{1}([-3,3]).$

As the PDE in (5.4) admits a closed-form solution parametric in $r$131313See
for example Section 4.5.1 of the lectures notes by Dr. Evy Kersalé
http://www1.maths.leeds.ac.uk/~kersale/2620/Notes/chapter_4.pdf., one can
easily bound $L_{1}(f)$, e.g., we use $L_{1}(f)=10$. Moreover, we set
$\sigma_{\xi}=10^{-12}$ to simulate numerical noise, set $\delta=10^{-6}$ and
perform the constrained optimization by means of Algorithm 4.1 (b) and by
using the potential function one can find for (5.4), that is, a function
$\varphi$ such that $u=\nabla\varphi$. Note that using our scheme and some
numerical PDE-solver as an inner-loop (instead of the closed-form solution) is
also possible, e.g., one needs to solve a linear system, not over
$\mathbb{R}$, but over $\mathbb{C}$. We select $8$ initial conditions
uniformly from $[1,8]$ and show the convergence in Figure 5.2141414See
http://wjongeneel.nl/PDE.gif for an animated version of Figure 2(c).. Note in
particular that the non-averaged iterates perform similar to their averaged
counterparts.

(a) Example 5.3, convergence of the cost (averaged estimator).

(b) Example 5.3, convergence of $r_{0}\in[1,2]$ (non-averaged).

(c) Example 5.3, typical result for $K=10^{6}$ (non-averaged).

Figure 5.2: PDE-constrained optimization. Each figure displays all available
data.

##  6 Discussion

###  6.1 On the necessity of leaving the real numbers

Given the results from the previous section, one might wonder if this
“complex-lifting” is needed. Real single-point gradient estimators evidently
exist, cf. [FKM04], but with problematic variance bounds for $\delta\downarrow
0$. The common solution is to bring back some relation with the (directional)
derivative [ADX10, NS17]. Hence, one might wonder if there is a purely real
analogue to (2.2). The next proposition strongly hints at a negative answer.

######  Proposition 6.1 (On the necessity of leaving the real numbers).

Consider some open, convex set $\mathcal{D}\subseteq\mathbb{R}^{n}$ with
$\mathrm{int}(\mathcal{D})\neq\emptyset$. Then, there does not exist a
continuous map $G:\mathbb{R}\to\mathbb{R}$ such that for all real-analytic
functions $f:\mathcal{D}\to\mathbb{R}$

$\lim_{\delta\downarrow 0}G\left(\tfrac{1}{\delta}f(x+\delta
y)\right)=Df(x)[y]=\langle\nabla f(x),y\rangle\quad\forall
x\in\mathcal{D},\,y\in\mathbb{S}^{n-1}.$

###### Proof.

As $f\in C^{\omega}(\mathcal{D})$ we can construct for sufficiently small
$\delta$ and any $y\in\mathbb{S}^{n-1}$ the convergent Taylor series of
$f(x+\delta y)$ around $x$ and as $G\in C^{0}(\mathbb{R})$ we can consider the
limit in $\delta$ with respect to the argument of $G$, hence, we have
$\lim_{\delta\downarrow 0}G\left(\tfrac{1}{\delta}f(x+\delta
y)\right)=G(\langle\nabla f(x),y\rangle+\lim_{\delta\downarrow
0}\tfrac{1}{\delta}f(x))$. When $f(x)=0$ we end up with the fixed-point
problem $G(\langle\nabla f(x),y\rangle)=\langle\nabla f(x),y\rangle$. As
$\mathcal{D}$ is open, convex and with
$\mathrm{int}(\mathcal{D})\neq\emptyset$, then for any $\alpha\in\mathbb{R}$
one can always find a pair $(f,y)$ such that $\langle\nabla
f(x),y\rangle=\alpha$, e.g., construct a linear function over $\mathcal{D}$.
Therefore, $G$ is forced to be the identity map on $\mathbb{R}$. Thereby,
obstructing the case $f(x)\neq 0$. ∎

Observe from the proof of Proposition 6.1 that if we would generalize
$G(\frac{1}{\delta}f(x+\delta y))$ to $\bar{G}(f(x+\delta y),\delta)$ with
$\bar{G}$ continuous in $\mathbb{R}\times\mathbb{R}_{>0}$, the conclusion
would not change.

###  6.2 $C^{\infty}$-smooth imaginary zeroth-order optimization

Consider the smooth function $\psi:\mathbb{R}\to\mathbb{R}$ defined by
$\psi(x)=|x|^{2}$. When evaluating $\psi$ at some complex point
$z=x+iy\in\mathbb{C}$ one finds that $\psi(z)=x^{2}+y^{2}$, as such, $\psi$
does not satisfy the Cauchy-Riemann equations and is nowhere (complex)
analytic. This, however, means that one cannot appeal to the complex-step
framework from [JYK21] cf. Section 3. Next, consider the prototypical smooth,
yet non-analytic, function $\varphi:\mathbb{R}\to\mathbb{R}$ defined by

$\varphi(x)=\begin{cases}\mathrm{exp}\left(\tfrac{-1}{x}\right)&\text{if
}x>0\\\ 0&\text{otherwise}\end{cases}.$

This function only fails to be analytic at $0$ and, interestingly, by the
following expansions of $\mathrm{exp}(-1/z)$

$\displaystyle\mathrm{exp}\left(\tfrac{-1}{x+iy}\right)=\mathrm{exp}\left(\tfrac{-x}{x^{2}+y^{2}}+i\tfrac{y}{x^{2}+y^{2}}\right)=\mathrm{exp}\left(\tfrac{-x}{x^{2}+y^{2}}\right)\left(\cos\left(\tfrac{y}{x^{2}+y^{2}}\right)+i\sin\left(\tfrac{y}{x^{2}+y^{2}}\right)\right)$
(6.1)

one can readily show that $\varphi$ does satisfy the Cauchy-Riemann equations.
Indeed, recall (2.2) and consider now the imaginary part of (6.1), then by the
series expansion of $\mathrm{exp}(\cdot)$ and $\sin(\cdot)$ one observes that

$\partial_{x}\,\mathrm{exp}\left(\tfrac{-1}{x}\right)=\tfrac{1}{x^{2}}\mathrm{exp}\left(\tfrac{-1}{x}\right)=\lim_{\delta\downarrow
0}\,\tfrac{1}{\delta}\mathrm{exp}\left(\tfrac{-x}{x^{2}+\delta^{2}}\right)\sin\left(\tfrac{\delta}{x^{2}+\delta^{2}}\right).$

Hence, although $\varphi\in C^{\infty}\setminus C^{\omega}$, the complex-step
framework is not obstructed.

It turns out that from a topological point of view, the function $\varphi$ is
somewhat of a special case. Let $\mathsf{X}$ be a topological space. Then the
set $M\subset\mathsf{X}$ is of the first category, in the sense of Baire, when
$M$ is a countable union of nowhere dense sets in $\mathsf{X}$. A set
$A\subseteq\mathsf{X}$ is said to be nowhere dense when $\mathrm{cl}(A)^{c}$
is dense in $\mathsf{X}$, or equivalently, when
$\mathrm{int}(\mathrm{cl}(A))=\emptyset$. Now one can show that under the sup-
norm, the complement to the space of nowhere differentiable functions in
$C^{0}([0,1])$ is of the first category [Fol99, Chapter 5]. Differently put,
almost every continuous function on $[0,1]$ is nowhere differentiable. A
similar topological statement can be made about nowhere analytic functions in
the space of smooth functions $C^{\infty}([0,1])$ under a sup-metric, e.g.,
see151515See in particular this post
https://web.archive.org/web/20161009194815/mathforum.org/kb/message.jspa?messageID=387148
by Dave L. Renfro for more context. [Dar73, Cat84]. Again, bluntly put, almost
every smooth function is nowhere analytic. An important question that comes
with such an observation is where in the space of smooth functions
optimization takes place?

###  6.3 Future work

This work exploits smoothness to be able to appeal to the Cauchy-Riemann
equations. Other work, like [PT90, BP16, APT20, NG21] exploit the knowledge of
smoothness and construct kernels to (optimally) filter out (all) low-order
errors. For increasing smoothness, however, we observe numerical instability
in this approach, that is, the kernels become ill-defined. It would be
worthwhile to further study how to exploit smoothness while taking the
implementation into consideration. Given Proposition 6.1, it would also be
interesting to explore the possibility of applying generalized versions of the
complex-step approach, e.g., using hyper-dual numbers to extract second-order
information [FA11].

This work is mostly positioned within the scope of randomized methods via
Lemma 3.3. Recent work indicated that in fact non-randomized methods can
outperform their randomized/smoothed counterparts [Ber+21, Sch22]. This
provides for interesting future work, especially in the presence of noise.
Estimating the noise statistics itself also provides for relevant future work
as it allows for a more appropriately scaled sequence of smoothing parameters.

###  6.4 Conclusion

We have presented a line of algorithms that can theoretically and practically
deal with any suitable sequence $\delta_{k}\to 0$ (conditioned on appropriate
stepsizes $\\{\mu_{k}\\}_{k}$). In contrast to [JYK21] we can also deal with
computational noise and demand less prior knowledge of problem parameters.

Only if we understand all the vulnerabilities of our algorithms — as esoteric
as they are — we can safely implement them. With that in mind, we hope this
work provides for more future work on numerical optimization.

## Appendix

This appendix contains auxiliary results related to the work above.

##  Appendix A Auxiliary results

The following results are well-known.

######  Lemma A.1 (Logarithm bound).

For any $J\in\mathbb{N}_{\geq 1}$ one has

$\textstyle\sum^{J}_{j=1}\tfrac{1}{j}\leq\log(J)+1.$ (A.1)

######  Lemma A.2 (Fractional bound).

For any $\beta\geq 1$ one has

$\textstyle\sum^{J}_{j=1}j^{-1+1/\beta}\leq\beta J^{1/\beta}.$ (A.2)

##  Appendix B Estimation of $\bar{\tau}(f)$ via diagonally dominant
programming

We highlight the basis pursuit approach as proposed in [AH17]. A constraint of
the form $P\in\mathsf{dd}^{n}$ translates to a set of linear constraints. The
same is true for $P_{z}\in\mathsf{dd}^{n}(U_{z})$ with

$\mathsf{dd}^{n}(U_{z})=\\{M\in\mathcal{S}^{n}:M=U_{z}^{\mathsf{T}}QU_{z},\,Q\in\mathsf{dd}^{n}\\}$

for some basis matrix $U_{z}$. Now to iteratively change the basis $U_{z}$ one
can use

$U_{z+1}=\mathsf{chol}(P_{z}),\quad U_{0}=I_{n},$

for $P_{z}$ the solution of the $z^{\mathrm{th}}$ program. By construction one
has $P_{z}\in\mathsf{dd}^{n}(U_{z+1})$ such that each new iteration is at
least as good as the previous one. As by the construction in Section 4.5 we
demand that $P_{z}\succeq\tau_{0}I_{n}\succ 0$, then, by [AH17, Theorem 3.1]
$P_{z}\to P^{\star}$ (weakly) for $z\to+\infty$ and $P^{\star}$ being the (a)
solution of the original problem. In practice, one could terminate the
algorithm when $Q_{z}$ is sufficiently close to $I_{n}$ and set
$\widehat{\tau}(f)=\min_{i}\\{\lambda_{i}(P_{z})\\}$, which can be found using
a dedicated large-scale algorithm.

To showcase the approach we redo Example 2.1, but by using (4.17). Here, we
fix a random pair $(A,b)$ and show for $100$ initial conditions
$x_{1}\overset{i.i.d.}{\sim}\mathcal{N}(0,I_{n})$ the effect of an improved
estimate of $\bar{\tau}$. Here, we apply the basis pursuit approach as
sketched above for $n$ iterations and set $\widehat{\tau}=\widehat{\tau}_{n}$.
The results are shown in Figure B.1. Again, we observe the benefit of
estimating $\bar{\tau}$, plus, we see that the inner-routine convergences
quickly, yet, usually from above. Quantifying the behaviour as seen in Figure
1(b) would be interesting and is left for future work.

(a) Algorithm 4.1 (b) in combination with DDP-based estimation of
$\bar{\tau}$.

(b) Basis pursuit convergence for the program (4.17).

Figure B.1: Numerical outcomes of the DDP example (Section B).

##  Appendix C Lipschitz inequalities

In this section we gather a variety of inequalities which come in useful
later. Note that convexity of $f$ is usually not a necessary assumption. If
$f$ is convex then, by [Nes03, Theorem 2.1.5] (1.2) implies that

$f(x)\geq f(y)+\langle\nabla f(y),x-y\rangle+\tfrac{1}{2L_{1}(f)}\|\nabla
f(x)-\nabla f(y)\|_{2}^{2},\quad\forall x,y\in\mathcal{D}$ (C.1)

and thus for any (local) minimum $x^{\star}$ such that $\nabla f(x^{\star})=0$
one has $2L_{1}(f)\left(f(x)-f(x^{\star})\right)\geq\|\nabla f(x)\|_{2}^{2}$.
Also, as [Nes11, Equation (6)], for $f\in C^{1,1}_{L_{1}(f)}(\mathcal{D})$ one
has

$|f(y)-f(x)-\langle\nabla
f(x),y-x\rangle|\leq\tfrac{1}{2}L_{1}(f)\|x-y\|_{2}^{2},\quad\forall
x,y\in\mathcal{D}.$ (C.2)

It follows from [Nes03, Lemma 1.2.4] that if $f\in C^{2}(\mathcal{D})$ then

$|f(y)-f(x)-\langle\nabla
f(x),y-x\rangle-\tfrac{1}{2}\langle\nabla^{2}f(x)(y-x),y-x\rangle|\leq\tfrac{1}{6}L_{2}(f)\|x-y\|_{2}^{3},\quad\forall
x,y\in\mathcal{D}.$ (C.3)

See that (1.3) is equivalent to

$|\langle\nabla^{2}f(x)u,u\rangle-\langle\nabla^{2}f(y)u,\rangle|\leq
L_{2}(f)\|x-y\|_{2}\quad\forall x,y\in\mathcal{D},\;u\in\mathbb{S}^{n-1},$
(C.4)

which is commonly referred to as $f$ being $3$rd-order smooth, cf. [BP16,
Section 1.1]. Now it follows directly from (C.3) and the definition of a
derivative that $f\in C^{3,2}_{L_{2}(f)}(\mathcal{D})$ implies that for all
$x\in\mathcal{D}$ one has

$\left|\partial^{3}_{t}f(x+tu)|_{t=0}\right|\leq L_{2}(f),\quad\forall
u\in\mathbb{S}^{n-1}.$ (C.5)

##  Appendix D Further numerical comments

Example 4.8 continued. In Figure 1(b) we see a clear difference in behaviour.
This can be explained by looking at the corresponding estimators. We see that
for the estimator as proposed in this work no cancellation occurs, while for
the frequently employed central-difference scheme as used in [APT20] the two
function evaluations can cancel catastrophically. See Figure 1(a)-1(b) and
Figure 1(c). We like to remark, in line with the analysis, that the scheme for
$\delta=1$ is better conditioned.

(a) CS ($\delta=1$) estimator

(b) CS ($\delta=10^{-100}$) estimator

(c) $\beta$ estimator

Figure D.1: Further remarks on Example 4.8.

All numerical experiments are carried out in MATLAB using the SDPT3 solver
[TTT99].

#### Data availability statement

All data generated or analysed during this study are included in this article.

#### Conflict of interest

The author has no competing interests to declare that are relevant to the
content of this article.

### Bibliography

## References

  * [Abr+18] Rafael Abreu, Zeming Su, Jochen Kamm and Jinghuai Gao  “On the accuracy of the Complex-Step-Finite-Difference method”  In _Journal of Computational and Applied Mathematics_ 340, 2018, pp. 390–403
  * [ADX10] Alekh Agarwal, Ofer Dekel and Lin Xiao  “Optimal algorithms for online convex optimization with multi-point bandit feedback.”  In _Conference on Learning Theory_ , 2010, pp. 28–40
  * [Aga+09] Alekh Agarwal, Martin J Wainwright, Peter Bartlett and Pradeep Ravikumar  “Information-theoretic lower bounds on the oracle complexity of convex optimization”  In _Neural Information Processing Systems_ , 2009, pp. 1–9
  * [AH17] Amir Ali Ahmadi and Georgina Hall  “Sum of squares basis pursuit with linear and second order cone programming”  In _Algebraic and Geometric Methods in Discrete Mathematics_ 685 American Mathematical Society Providence, RI, 2017, pp. 27–53
  * [AH17a] Charles Audet and Warren Hare  “Derivative-free and Blackbox Optimization”  Springer, 2017
  * [AL81] Helmut Alt and Jan Leeuwen  “The complexity of basic complex operations”  In _Computing_ 27.3 Springer, 1981, pp. 205–215
  * [AM19] Amir Ali Ahmadi and Anirudha Majumdar  “DSOS and SDSOS optimization: more tractable alternatives to sum of squares and semidefinite optimization”  In _SIAM Journal on Applied Algebra and Geometry_ 3.2 SIAM, 2019, pp. 193–230
  * [AMA05] Pierre-Antoine Absil, Robert Mahony and Benjamin Andrews  “Convergence of the iterates of descent methods for analytic cost functions”  In _SIAM Journal on Optimization_ 16.2 SIAM, 2005, pp. 531–547
  * [AMH10] Awad Al-Mohy and Nicholas Higham  “The complex step approximation to the Fréchet derivative of a matrix function”  In _Numerical Algorithms_ 53, 2010, pp. 133–148
  * [APT20] Arya Akhavan, Massimiliano Pontil and Alexandre Tsybakov  “Exploiting higher order smoothness in derivative-free optimization and continuous bandits”  In _Neural Information Processing Systems_ , 2020, pp. 9017–9027
  * [AS21] Ahmad Ajalloeian and Sebastian U Stich  “On the convergence of SGD with biased gradients”, 2021  arXiv:2008.00051
  * [ASM15] Rafael Abreu, Daniel Stich and Jose Morales  “The Complex-Step-Finite-Difference method”  In _Geophysical Journal International_ 202.1, 2015, pp. 72–93
  * [Ber+21] Albert S Berahas, Liyuan Cao, Krzysztof Choromanski and Katya Scheinberg  “A theoretical and empirical comparison of gradient approximations in derivative-free optimization”  In _Foundations of Computational Mathematics_ Springer, 2021, pp. 1–54
  * [BG21] Krishnakumar Balasubramanian and Saeed Ghadimi  “Zeroth-order nonconvex stochastic optimization: handling constraints, high dimensionality, and saddle points”  In _Foundations of Computational Mathematics_ Springer, 2021, pp. 1–42
  * [BP16] Francis Bach and Vianney Perchet  “Highly-smooth zero-th order online optimization”  In _Conference on Learning Theory_ , 2016, pp. 257–283
  * [Cat84] FS Cater  “Differentiable, nowhere analytic functions”  In _The American Mathematical Monthly_ 91.10, 1984, pp. 618–624
  * [CH04] MG Cox and PM Harris  “Software support for meteorology best practice guide no. 11”, 2004
  * [Che88] Hung Chen  “Lower rate of convergence for locating a maximum of a function”  In _The Annals of Statistics_ 16.3, 1988, pp. 1330–1334
  * [CSV09] Andrew R. Conn, Katya Scheinberg and Luis N. Vicente  “Introduction to Derivative-Free Optimization”  Society for IndustrialApplied Mathematics, 2009
  * [CTL22] Xin Chen, Yujie Tang and Na Li  “Improve single-point zeroth-order optimization using high-pass and low-pass filters”  In _International Conference on Machine Learning_ , 2022, pp. 3603–3620
  * [CWF20] Charles Champagne Cossette, Alex Walsh and James Richard Forbes  “The complex-step derivative approximation on matrix lie groups”  In _IEEE Robotics and Automation Letters_ 5.2, 2020, pp. 906–913
  * [d’A08] Alexandre d’Aspremont  “Smooth optimization with approximate gradient”  In _SIAM Journal on Optimization_ 19.3 SIAM, 2008, pp. 1171–1183
  * [Dar73] RB Darst  “Most infinitely differentiable functions are nowhere analytic”  In _Canadian Mathematical Bulletin_ 16.4, 1973, pp. 597–598
  * [DGN14] Olivier Devolder, François Glineur and Yurii Nesterov  “First-order methods of smooth convex optimization with inexact oracle”  In _Mathematical Programming_ 146.1 Springer, 2014, pp. 37–75
  * [Duc+15] John C Duchi, Michael I Jordan, Martin J Wainwright and Andre Wibisono  “Optimal rates for zero-order convex optimization: The power of two function evaluations”  In _IEEE Transactions on Information Theory_ 61.5 IEEE, 2015, pp. 2788–2806
  * [FA11] Jeffrey Fike and Juan Alonso  “The development of Hyper-Dual numbers for exact second-derivative calculations”  In _49th AIAA Aerospace Sciences Meeting_ , 2011
  * [Fab71] V Fabian  “Stochastic approximation”  In _Symposium on Optimizing Methods in Statistics (1971: Ohio State University)_ , 1971, pp. 439–470  Academic Press
  * [Faz+18] Maryam Fazel, Rong Ge, Sham Kakade and Mehran Mesbahi  “Global convergence of policy gradient methods for the Linear Quadratic regulator”  In _International Conference on Machine Learning_ , 2018, pp. 1467–1476
  * [FKM04] Abraham Flaxman, Adam Tauman Kalai and H. Brendan McMahan  “Online convex optimization in the bandit setting: gradient descent without a gradient”, 2004  arXiv:0408007
  * [Fol99] Gerald B. Folland  “Real Analysis”  Wiley-Interscience, 1999
  * [Gas+17] Alexander V Gasnikov, Ekaterina A Krymova, Anastasia A Lagunovskaya, Ilnura N Usmanova and Fedor A Fedorenko  “Stochastic online optimization. Single-point and multi-point non-linear multi-armed bandits. Convex and strongly-convex case”  In _Automation and remote control_ 78.2 Springer, 2017, pp. 224–234
  * [GHL04] Sylvestre Gallot, Dominique Hulin and Jacques Lafontaine  “Riemannian Geometry”  Springer, 2004
  * [GL13] Saeed Ghadimi and Guanghui Lan  “Stochastic first-and zeroth-order methods for nonconvex stochastic programming”  In _SIAM Journal on Optimization_ 23.4 SIAM, 2013, pp. 2341–2368
  * [HL14] Elad Hazan and Kfir Levy  “Bandit convex optimization: towards tight bounds”  In _Neural Information Processing Systems_ , 2014, pp. 784–792
  * [HRB08] Elad Hazan, Alexander Rakhlin and Peter Bartlett  “Adaptive Online Gradient Descent”  In _Advances in Neural Information Processing Systems_ , 2008, pp. 65–72
  * [HS21] Warren Hare and Kashvi Srivastava  “A numerical study of applying complex-step gradient and Hessian approximations in blackbox optimization”, 2021  URL: https://www.researchgate.net/publication/355081274_A_Numerical_Study_of_Applying_Complex-step_Gradient_and_Hessian_Approximations_in_Blackbox_Optimization
  * [Hu+16] Xiaowei Hu, LA Prashanth, András György and Csaba Szepesvari  “(Bandit) convex optimization with biased noisy gradient oracles”  In _Artificial Intelligence and Statistics_ , 2016, pp. 819–828
  * [JNR12] Kevin G Jamieson, Robert Nowak and Ben Recht  “Query complexity of derivative-free optimization”  In _Neural Information Processing Systems_ , 2012, pp. 2681–2689
  * [JYK21] Wouter Jongeneel, Man-Chung Yue and Daniel Kuhn  “Small errors in random zeroth-order optimization are imaginary”, 2021  arXiv:2103.05478
  * [KC78] Harold Joseph Kushner and Dean S Clark  “Stochastic Approximation Methods for Constrained and Unconstrained Systems”  Springer, 1978
  * [Kra00] Steven G. Krantz  “Function Theory of Several Complex Variables”  AMS Chelsea Publishing, 2000
  * [KSST09] Sham M Kakade, Shai Shalev-Shwartz and Ambuj Tewari  “On the duality of strong convexity and strong smoothness: Learning applications and matrix regularization”, 2009  URL: https://home.ttic.edu/~shai/papers/KakadeShalevTewari09.pdf
  * [KW52] Jack Kiefer and Jacob Wolfowitz  “Stochastic estimation of the maximum of a regression function”  In _The Annals of Mathematical Statistics_ , 1952, pp. 462–466
  * [KY03] Harold Kushner and G George Yin  “Stochastic approximation and recursive algorithms and applications”  Springer Science & Business Media, 2003
  * [Lee13] John M. Lee  “Introduction to Smooth Manifolds”  Springer, 2013
  * [Liu+20] S. Liu, P. Y. Chen, B. Kailkhura, G. Zhang, A. O. Hero III and P. K. Varshney  “A primer on zeroth-order optimization in signal processing and machine learning: principals, recent advances, and applications”  In _IEEE Signal Processing Magazine_ 37.5, 2020, pp. 43–54
  * [LLZ21] Henry Lam, Haidong Li and Xuhui Zhang  “Minimax efficient finite-difference stochastic gradient estimators using black-box function evaluations”  In _Operations Research Letters_ 49.1 Elsevier, 2021, pp. 40–47
  * [LM67] J. N. Lyness and C. B. Moler  “Numerical differentiation of analytic functions”  In _SIAM Journal on Numerical Analysis_ 4.2, 1967, pp. 202–210
  * [LMW19] Jeffrey Larson, Matt Menickelly and Stefan M. Wild  “Derivative-free optimization methods”  In _Acta Numerica_ 28 Cambridge University Press, 2019, pp. 287–404
  * [LT93] Zhi-Quan Luo and Paul Tseng  “Error bounds and convergence analysis of feasible descent methods: a general approach”  In _Annals of Operations Research_ 46.1 Springer, 1993, pp. 157–178
  * [Mal+19] Dhruv Malik, Ashwin Pananjady, Kush Bhatia, Koulik Khamaru, Peter Bartlett and Martin Wainwright  “Derivative-free methods for policy optimization: Guarantees for linear quadratic systems”  In _International Conference on Artificial Intelligence and Statistics_ , 2019, pp. 2916–2925
  * [MHA20] Anirudha Majumdar, Georgina Hall and Amir Ali Ahmadi  “Recent scalability improvements for semidefinite programming with applications in machine learning, control, and robotics”  In _Annual Review of Control, Robotics, and Autonomous Systems_ 3, 2020, pp. 331–360
  * [Moc12] Jonas Mockus  “Bayesian Approach to Global Optimization: Theory and Applications”  Springer Science & Business Media, 2012
  * [Mor58] Charles B Morrey  “On the analyticity of the solutions of analytic non-linear elliptic systems of partial differential equations: Part I. Analyticity in the interior”  In _American Journal of Mathematics_ 80.1, 1958, pp. 198–218
  * [Mor58a] Charles B Morrey  “On the analyticity of the solutions of analytic non-linear elliptic systems of partial differential equations: Part II. Analyticity at the boundary”  In _American Journal of Mathematics_ 80.1, 1958, pp. 219–237
  * [MSA03] Joaquim R. R. A. Martins, Peter Sturdza and Juan J. Alonso  “The Complex-Step derivative approximation”  In _ACM Trans. Math. Softw._ 29.3, 2003, pp. 245––262
  * [Nes03] Yurii Nesterov  “Introductory Lectures on Convex Optimization: a Basic Course”  Springer Science & Business Media, 2003
  * [Nes11] Yurii Nesterov  “Random gradient-free minimization of convex functions”, 2011
  * [Nes13] Yu Nesterov  “Gradient methods for minimizing composite functions”  In _Mathematical Programming_ 140.1 Springer, 2013, pp. 125–161
  * [NG21] Vasilii Novitskii and Alexander Gasnikov  “Improved exploiting higher order smoothness in derivative-free optimization and continuous bandit”, 2021  arXiv:2101.03821
  * [NS17] Yurii Nesterov and Vladimir Spokoiny  “Random gradient-free minimization of convex functions”  In _Foundations of Computational Mathematics_ 17.2 Springer, 2017, pp. 527–566
  * [NS18] Filip Nikolovski and Irena Stojkovska  “Complex-step derivative approximation in noisy environment”  In _Journal of Computational and Applied Mathematics_ 327, 2018, pp. 64–78
  * [NY83] Arkadi Semenovich Nemirovsky and David Borisovich Yudin  “Problem Complexity and Method Efficiency in Optimization”  Wiley, 1983
  * [OC15] Brendan O’donoghue and Emmanuel Candes  “Adaptive restart for accelerated gradient schemes”  In _Foundations of computational mathematics_ 15.3 Springer, 2015, pp. 715–732
  * [Oss78] Robert Osserman  “The isoperimetric inequality”  In _Bulletin of the American Mathematical Society_ 84.6, 1978, pp. 1182–1238
  * [Ove01] Michael L Overton  “Numerical Computing with IEEE Floating Point Arithmetic”  Society for IndustrialApplied Mathematics, 2001
  * [Pol86] J.W. Polderman  “A note on the structure of two subsets of the parameter space in adaptive control problems”  In _Systems & Control Letters_, 1986, pp. 25–34
  * [PS17] Boris Polyak and Pavel Shcherbakov  “Why does Monte Carlo fail to work properly in high-dimensional optimization problems?”  In _Journal of Optimization Theory and Applications_ 173.2 Springer, 2017, pp. 612–627
  * [PT90] Boris Teodorovich Polyak and Aleksandr Borisovich Tsybakov  “Optimal order of accuracy of search algorithms in stochastic optimization”  In _Problemy Peredachi Informatsii_ 26.2, 1990, pp. 45–53
  * [Riv07] Igor Rivin  “Surface area and other measures of ellipsoids”  In _Advances in Applied Mathematics_ 39.4 Elsevier, 2007, pp. 409–427
  * [RSS12] Alexander Rakhlin, Ohad Shamir and Karthik Sridharan  “Making gradient descent optimal for strongly convex stochastic optimization”  In _International Conference on Machine Learning_ , 2012, pp. 1571–1578
  * [Sch22] Katya Scheinberg  “Finite Difference Gradient Approximation: To Randomize or Not?”  In _INFORMS Journal on Computing_ INFORMS, 2022, pp. 1–5
  * [Sha13] Ohad Shamir  “On the complexity of bandit and derivative-free stochastic convex optimization”  In _Conference on Learning Theory_ , 2013, pp. 3–24
  * [Sha17] Ohad Shamir  “An optimal algorithm for bandit and zero-order convex optimization with two-point feedback”  In _The Journal of Machine Learning Research_ 18.1, 2017, pp. 1703–1713
  * [Shi+21] Hao-Jun Michael Shi, Melody Qiming Xuan, Figen Oztoprak and Jorge Nocedal  “On the numerical performance of derivative-free optimization methods based on finite-difference approximations”, 2021  arXiv:2102.09762
  * [Spa05] James C Spall  “Introduction to Stochastic Search and Optimization: Estimation, Simulation, and Control”  John Wiley & Sons, 2005
  * [ST98] William Squire and George Trapp  “Using complex variables to estimate derivatives of real functions”  In _SIAM Review_ 40.1, 1998, pp. 110–112
  * [TSAK21] Bahar Taşkesen, Soroosh Shafieezadeh-Abadeh and Daniel Kuhn  “Semi-discrete optimal transport: hardness, regularization and numerical solution”, 2021  arXiv:2103.06263
  * [TTT99] Kim-Chuan Toh, Michael J Todd and Reha H Tütüncü  “SDPT3—a MATLAB software package for semidefinite programming, version 1.3”  In _Optimization methods and software_ 11.1-4 Taylor & Francis, 1999, pp. 545–581
  * [WS21] Long Wang and James C Spall  “Improved SPSA using complex variables with applications in optimal control problems”  In _American Control Conference_ , 2021, pp. 3519–3524  IEEE
  * [WZS21] Long Wang, Jingyi Zhu and James C Spall  “Model-free optimal control using SPSA with complex variables”  In _Annual Conference on Information Sciences and Systems_ , 2021, pp. 1–5  IEEE
  * [Zha+22] Yan Zhang, Yi Zhou, Kaiyi Ji and Michael M. Zavlanos  “A new one-point residual-feedback oracle for black-box learning and control”  In _Automatica_ 136.C, 2022