# Small Errors in Random Zeroth Order Optimization
are Imaginary

Wouter Jongeneela Man-Chung Yueb Daniel Kuhna

( aRisk Analytics and Optimization Chair, École Polytechnique Fédérale de
Lausanne<EMAIL_ADDRESS>
bDepartment of Applied Mathematics, The Hong Kong Polytechnic University,
<EMAIL_ADDRESS>

)

###### Abstract

The vast majority of zeroth order optimization methods try to imitate first
order methods via some smooth approximation of the gradient. Here, the smaller
the smoothing parameter, the smaller the gradient approximation error. We show
that for the majority of zeroth order methods this smoothing parameter can
however not be chosen arbitrarily small as numerical cancellation errors will
dominate. As such, theoretical and numerical performance could differ
significantly. Using classical tools from numerical differentiation we will
propose a new smoothed approximation of the gradient that can be integrated
into general zeroth order algorithmic frameworks. Since the proposed smoothed
approximation does not suffer from cancellation error, the smoothing parameter
(and hence the approximation error) can be made arbitrarily small. Sublinear
convergence rates for algorithms based on our smoothed approximation are
proved. Numerical experiments are also presented to demonstrate the
superiority of algorithms based on the proposed approximation.

##  1 Introduction

Let $f:\mathcal{D}\subseteq\mathbb{R}^{n}\to\mathbb{R}$ be some objective
function, $\mathcal{D}$ an open set and $\mathcal{X}\subseteq\mathcal{D}$ a
closed non-empty set, then, we are interested in solving problems of the form

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

Assuming some optimizer in (1.1) exists, we denote any by $x^{\star}$. In
solving (1.1) we will not have a standard arsenal of tools at our disposal,
instead, we assume to have access to function evaluations only, that is, a
finite sequence of the form $f(x_{0}),f(x_{1}),\dots,f(x_{K})$. In other
words, the objective function is unknown to us, but sampling its actions is
possible. This approach falls under the umbrella of zeroth order optimization,
or equivalently derivative-free optimization.

In motivating this setting we briefly highlight some remarks taken from the
monograph [CSV09] and the paper [Nes11]. Zeroth order optimization is a
framework intended to deal with problems (1.1) for which first order methods
do not (easily) apply. For example, obtaining a reasonable approximation of
the gradient $\nabla f(x)$ might be costly or simply impossible. The latter
scenario is true for many simulation-based problems where a closed-form
expression of $f$ is unavailable, e.g., for many modern engineering problems
$x$ relates to a design parameter and $f$ to an online physics simulator only
capable of providing the evaluation $f(x)$. An equally complicated setting
could be that of $f$ being itself of a variational form, e.g., in minimax
optimization problems $f$ is of the form $f(x)=\sup_{z\in\mathcal{Z}}h(x,z)$.
Explicit gradients are usually not available, or worse, the inner maximization
can only be solved numerically. Another setting frequently encountered in
practice is that of finding the optimal values of tuning parameters in
numerical algorithms. The dependence of the solution on this parameter is
almost always implicit which hinders classical gradient descent. However, as
function evaluations are possible, zeroth order methods are possible.
Furthermore, as only function evaluations are used, this approach is rather
crude and it is believed that this property makes it suitable for non-convex
problems, e.g., to escape saddle-points.

There are essentially two lines of attack when it comes to zeroth order
optimization, one could use the available data to construct an explicit
$C^{r}$ local model of $f$ for $r\geq 1$, say $\tilde{f}$ and apply
appropriate $r^{\mathrm{th}}$ order optimization tools to locally improve the
current iterate. Clearly, if $\tilde{f}$ is of high quality, then this is a
rather powerful method. However, to construct such a $\tilde{f}$ a lot of data
is usually needed.

To address this, a different and increasingly popular approach — which is also
the topic of this paper — is to approximate the gradient of $f$ directly from
the data. To do so, the vast majority of these zeroth order optimization
methods introduce a $\delta$-smoothed version of $f$, denoted $f_{\delta}$,
for $f_{\delta}$ being $\delta$-close to $f$. The idea is to construct an
unbiased estimator $g_{\delta}(x)$ for $\nabla f_{\delta}(x)$ via sampling the
objective $f$ at selected data $\\{x_{k}\\}_{k=0}^{K}$ around $x$. As one can
usually show that $\|\nabla f_{\delta}(x)-\nabla f(x)\|\leq O(\delta^{p})$ for
some $p\geq 1$, the estimator $g_{\delta}(x)$ can be used as a surrogate
gradient in for example a gradient-descent algorithm. Compared to model-based
methods, the key benefit is that less data and less assumptions are needed to
find an approximation of $x^{\star}$ in (1.1). Nevertheless, this usually
comes at the cost of weaker guarantees. Clearly, one likes to select the
smallest possible smoothing parameter $\delta$ as this increases approximation
quality and thereby presumably performance of the algorithm. Unfortunately,
the numerically useful range of $\delta$ — and thereby of the full scheme — is
usually unknown or at least not discussed.

In this work we will exploit insights from numerical differentiation to
propose a new smoothed approximation $f_{\delta}$ such that $\delta$, and
thereby the error in the gradient approximation scheme, can be made
arbitrarily small, which is not possible using conventional schemes.
Algorithms based on our new smoothed approximation enjoy provable convergence
rates, are more numerically stable, and most importantly, outperform those
algorithms based on existing smoothed approximation in terms of both accuracy
and computational speed as shown in our experiments.

#### Notation

We denote the real and imaginary parts of a complex number $z=a+ib$ by
$\Re(z)=a$ and $\Im(z)=b$, respectively. Let
$\mathbb{B}^{n}:=\\{x\in\mathbb{R}^{n}:\|x\|_{2}\leq 1\\}$ and
$\mathbb{S}^{n-1}:=\\{x\in\mathbb{R}^{n}:\|x\|_{2}=1\\}$ be the unit Euclidean
ball and sphere, respectively, in $\mathbb{R}^{n}$. Throughout the paper, if
distributions are declared to be uniform over $\mathcal{Y}$, we simply write
$y\sim\mathcal{Y}$ and additionally

$\mathbb{E}_{y\sim\mathcal{Y}}[g(y)]=\frac{1}{\mathrm{vol}(\mathcal{Y})}\int_{\mathcal{Y}}g(y)dy,$

with $dy$ denoting the Lebesgue measure of $\mathcal{Y}$. If
$f:\mathcal{D}\to\mathbb{R}$ is $r$ times continuously differentiable, we
denote this by $f\in C^{r}(\mathcal{D})$. Given two real-valued funtions
$h_{1}$ and $h_{2}$, then if there exists a function $h_{3}$ such that
$h_{1}(x)\leq h_{2}(x)+h_{3}(x)$ with $\lim_{x\to\infty}h_{3}(x)/h_{1}(x)=0$
we denote this asymptotic inequality by $h_{1}\lesssim h_{2}$. At last,
$C_{1},C_{2},\dots,C_{y},C_{z}$ always denote non-negative universal constants
which can change from line to line.

###  1.1 Related work

Zeroth order optimization based on a single-point oracle was arguably for the
first time described in [NY83, Section 9.3]. There, the authors approximate a
function $f:\mathbb{R}^{n}\to\mathbb{R}$ by the smoothed objective

$f_{\delta}(x)=\frac{1}{\mathrm{vol}(\mathbb{B}^{n})}\int_{\mathbb{B}^{n}}f(x+\delta
v)dv,$ (1.2)

where the amount of smoothing is controlled by the smoothing parameter
$\delta>0$. Then, the key observation111See the proof of Lemma 3.4 for more on
how to pass from (1.2) to (1.3). is that $\nabla f_{\delta}(x)$ admits an
integral representation of $f$ as well

$\nabla
f_{\delta}(x)=\frac{n}{\delta}\frac{1}{\mathrm{vol}(\mathbb{S}^{n-1})}\int_{\mathbb{S}^{n-1}}f(x+\delta
u)udu.$ (1.3)

Hence, a natural one-point candidate to approximate the gradient of $f$ would
be

$g_{\delta}(x)=\frac{n}{\delta}f(x+\delta u)u,\quad u\sim\mathbb{S}^{n-1}$
(1.4)

as $\mathbb{E}_{u\sim\mathbb{S}^{n-1}}[g_{\delta}(x)]=\nabla f_{\delta}(x)$.
The authors in [FKM04] applied this approach in the bandit context. However,
already in [NY83] it was remarked that the variance of (1.4) is problematic
for $\delta\downarrow 0$, which is inconvenient as a smaller $\delta$ in (1.2)
leads to a smaller bias. To overcome this deficiency, the key observation from
[ADX10, Nes11] is then to give (1.4) again the interpretation of a directional
derivative by stressing that one should consider a multi-point oracle of the
form

$g^{\prime}_{\delta}(x)=\frac{n}{\delta}\left(f(x+\delta u)-f(x)\right)u,\quad
u\sim\mathbb{S}^{n-1}.$ (1.5)

Although the expected values of $g_{\delta}(x)$ and $g^{\prime}_{\delta}(x)$
are equivalent, their variance is different since (1.5) is bounded in $\delta$
while (1.4) is not. One problem pertains, $\delta$ needs to be made
arbitrarily small to solve the minimization problem up to arbitrary precision,
yet for sufficiently small $\delta$, $f(x+\delta u)$ and $f(x)$ will be
numerically indistinguishable and hence oracles of the form (1.5) are still
prone to numerical complications and cannot always provide arbitrarily high
precision indeed. When function evaluations are noisy this is even more true
[Lia+16]. Also, see [HL14] for further generalizations of the oracle (1.4) and
[Duc+15, LLZ21] for more on the optimality of multi-point schemes.

We assume that $f$ is fixed and that the function measurements are
deterministic. Using deterministic function evaluations one could try to
directly approximate $\nabla f(x)$, e.g., by a coordinate-wise finite-
difference method. This line of work relates to inexact gradient methods
[d’A08, DGN14]. However, as in [Nes11, NS17], we will use a randomized
(analysis) method. One of the benefits is that we have a clear framework to
understand convergence of our algorithm, namely, we will consider convergence
in expectation $f(x_{k})\to f(x^{\star})$ whereas in fully deterministic
inexact methods a bias inevitably prevails.

Although we use a randomized method, as in [NS17] we assume to have access to
a deterministic oracle. This is in contrast with a variety of other works
which assume to have access only to noisy function evaluations of the form
$f(x)+\xi$, $\xi\overset{i.i.d.}{\sim}\Xi$. We note however that most methods
assume the noise to be well-behaved in the sense that its distribution is
light-tailed and hence the uncertainties can be well handled by averaging. If
the function evaluation contains adversarial noise, negative results are
obtained in [SV15].

Although our function measurements are without noise, we will randomly sample
function evaluations and hence we are interested in construction a sequence of
iterates $x_{0},x_{1},\dots,x_{K-1}$ which minimize the following regret

$\frac{1}{K}\sum^{K-1}_{k=0}\mathbb{E}\left[f(x_{k})-f(x^{\star})\right].$
(1.6)

Then, if the regret (1.6) is bounded by some $\epsilon\geq 0$ and $f$ is
convex, the averaged estimator $\bar{x}_{K-1}:=(1/K)\sum^{K-1}_{k=0}x_{k}$
achieves the expected suboptimality gap (optimization error)

$\mathbb{E}\left[f(\bar{x}_{K-1})-f(x^{\star})\right]\leq\frac{1}{K}\sum^{K-1}_{k=0}\mathbb{E}\left[f(x_{k})-f(x^{\star})\right]\leq\epsilon.$
(1.7)

We are interested in this suboptimality gap $\epsilon$ as a function of the
samplesize $K$ and a variety of other problem parameters such as $n$, the
dimension of $x$, and the Lipschitz constants of the derivatives of $f$. If
the problem setting allows for it, we analyze the gap
$\mathbb{E}[f(x_{k})-f(x^{\star})]$ as well, e.g., when $f$ is strongly
convex. If the samplesize $K$ appears in the suboptimality gap one usually
talks about the convergence rate, e.g., $\epsilon=O(1/K)$.

In [Nes11, NS17] the smoothing is carried out via passing to a Gaussian kernel
in combination with a multi-point oracle of the form (1.5). The authors
conclude with stating that random methods are typically $O(n)$ times slower
than their deterministic counterparts. For example, if $f$ is convex and has a
Lipschitz-continuous gradient, then the rate provided by [NS17] is $O(n/K)$.
In [Duc+15, Sha17] these results are improved, conditioned on having access to
multiple function evaluations. See also [Gas+17]. [GL13] consider extending
the nonconvex programming part of [Nes11] and [BG18] look at the method
proposed by [Nes11] via the Gaussian Stein identity and show that
approximating second order information is also possible.

Simultaneously, in [Sha13], the author provides lower bounds and shows that in
the derivative-free case, most errors cannot vanish faster than
$O(n/\sqrt{K})$. However, under certain conditions like strong-convexity, a
$O(n/K)$ rate is possible. Recall however that in [BM13] it is shown that by
exploiting problem structure $O(1/K)$ rates are possible for some non-strongly
convex problems.

Although our focus is not on the stochastic oracle setting we highlight
another line of closely related work. Starting with [PT90], a variety of
authors looked into exploiting higher-order smoothness. For example, in [BP16]
the ideas from [PT90] are extended to the online case and in [APT20, NG21] the
convergence rates are subsequently improved under increasingly weak noise
assumptions. It is imperative to remark that in the stochastic setting the
selection of the smoothing parameter $\delta$ has a twofold motivation. On the
one hand the estimation error should be small, but on the other hand a small
variance is preferred and usually this variance has a term proportional to
$1/\delta$. An increasingly popular method to further reduce the variance is
to make use of mini-batch methods, e.g., see [Ji+19].

Also, recent work considers extending zeroth order optimization ideas to the
Riemannian setting [LBM20] and even a more limited supply of information can
be considered [Cai+21].

One of the appealing aspects of zeroth order optimization is its simplicity.
As such we present simple algorithms with pre-determined stepsizes, but note
that line-search (direct search) and model-based methods provide for an
interesting alternative outlook [CSV09, SMG13, BCS19], e.g., to exploit
underlying low-dimensional structures [Gol+19].

For more information see the recent surveys on zeroth order (derivative-free)
optimization [LMW19, Liu+20].

###  1.2 Numerical considerations

The majority of the zeroth order optimization literature is approximating a
variety of derivatives and indeed gradients but unfortunately usually without
taking numerical considerations into account. For example, given some $f\in
C^{r}(\mathbb{R})$ with $r\geq 1$, then in the approximation

$\partial_{x}f(x)\approx\frac{f(x+\delta)-f(x)}{\delta}$ (1.8)

one cannot make $\delta>0$ arbitrarily small and expect to recover
$\partial_{x}f(x)$. Theoretically this is correct, but on a machine,
$f(x+\delta)$ and $f(x)$ coincide for sufficiently small $\delta>0$. For
example, let $f:\mathbb{R}\to\mathbb{R}$ be locally Lipschitz around $x$, then
there is some constant $L$ such that $|f(x+\delta)-f(x)|\leq L\cdot|\delta|$
for sufficiently small $\delta$. Now we see that if $L\cdot|\delta|$ is
smaller than machine-precision, numerical cancellation errors are prone to
dominate, see [Ove01, Chapter 11]. Running into these machine-precision
problems is inherent to finite-difference (multi-point) methods like (1.8),
yet, the vast majority of zeroth order oracles is of this form whilst
demanding $\delta$ to be arbitrarily small, e.g., in [APT20, Theorem 3.1] the
authors require $\delta_{k}=O(1/\sqrt{k})$ on a multi-point scheme. In a
small-noise regime this will be eventually problematic.

The pioneering and celebrated construction as proposed by [NY83] alleviates —
to the best of our knowledge unintentionally — part of this numerical issue by
approximating the gradient of $f$ using an oracle of the form
$\frac{n}{\delta}f(x+\delta u)u$, $u\sim\mathbb{S}^{n-1}$. Since there is just
one function call (single-point), the oracle has the potential to be
numerically more stable. However, as explained before, the variance is ill-
behaved. Therefore, we are motivated to propose a numerically stable oracle,
yet, allowing for convergence rates comparable to the modern literature. We
will propose a special one-point oracle which inherits the numerical
properties of one-point oracles, i.e., no cancellation errors, at the same
time, this oracle will display the approximation quality and convergence speed
of two-point oracles. All of this is achieved without introducing significant
computational overhead.

Throughout we make two assumptions which deserve to be emphasized.

######  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})$.

Although somewhat non-standard in the optimization literature, looking for
example at zeroth order optimization in the context of reinforcement learning
[Faz+18, Mal+19], we like to remark that for these particular policy iteration
problem the objective is real-analytic in the feedback gain [Pol86]. Also note
that $f\in C^{\omega}(\mathbb{R}^{n})$ does not immediately imply
$\beta^{\mathrm{th}}$ order smoothness as used in for example [BP16], e.g.,
$0^{\mathrm{th}}$ order smoothness translates to $f$ being bounded and
$1^{\mathrm{st}}$ order smoothness translates to $f$ being (globally)
Lipschitz.

Assumption 1.1 is mainly there to allow for the next assumption.

######  Assumption 1.2 (Complex oracle).

Consider some unknown function $f\in C^{\omega}(\mathcal{D})$, we assume to
have access to an oracle which can output $f(x+iy)\in\mathbb{C}$, for any
$x+iy=:z\in\Omega\subseteq\mathbb{C}^{n}$ with $\mathcal{D}\subset\Omega$.

Hence, the oracle must be able to process complex computations, plus, for
simplicity — and to some extent realism — we assume it to be deterministic.
The set $\Omega$ will be specified later on.

#### Contribution

It was pointed out in [Duc+15, LLZ21] that, when looking at theoretical
convergence rates, multi-point schemes are preferred in zeroth order
optimization over algorithms using single function evaluations. However, from
a numerical point of view, these single-point schemes are inherently more
robust as in contrast to multi-point schemes cancellation errors can be
avoided. In this work we combine the fast convergence rate and low variance of
multi-point schemes with the numerical robustness of the single-point
approach. All of this is made possible via passing to the complex domain.

Let $R:=\|x_{0}-x^{\star}\|_{2}$, $F:=f(x_{0})-f(x^{\star})$, then, under the
assumption that $f$ is real-analytic and has a $L$-Lipschitz continuous
gradient, we show for the convex case (unconstrained and constrained) that the
optimization error scales as $O(nLR^{2}/K)$. For the $\tau$-strongly convex
case, the optimization error decays at a global linear rate of the form
$O((1-\tau/(4nL))^{K}LR^{2})$, whereas for the non-convex case a local
$O(nLF/K)$ rate is shown. All these results are similar to the rates as in
[NS17]. In fact, these rates are sharper than the noise-free rates as provided
in [APT20], where higher-order smoothness — although for the perturbed setting
— is explicitly taken into account. The key difference with all of the
aforementioned works is that we can select our smoothing parameter as being
for example of the form $\delta_{k}=\delta/(1+k)$ without the fear of
numerical problems.

What is more, 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 a multi-
point method. This observation motivates [Zha+20] to use some form of memory
such that their oracle only demands a single new point each call.
Nevertheless, their oracle results in a multi-point method. Although their
algorithm is faster than [FKM04], it is slower than [NS17]. On the contrary,
we will provide the first real single-point method which is as fast, but
usually faster, than the multi-point methods by [NS17]. Numerical experiments
will corroborate this claim.

#### Structure

In Section 2 we highlight the complex-step method from numerical
differentiation along with indispensable tools from Complex analysis. Then, in
Section 3 we introduce our complex-step oracle and in Sections 4-6 we analyze
the performance of this oracle in the context of convex, strongly convex and
non-convex optimization algorithms, respectively. Then, in Section 7 we
provide a variety of numerical experiments to substantiate the claims. At
last, a few auxiliary results are relegated to the Appendix.

##  2 Preliminaries

As we assume our oracle to handle complex computations we first highlight what
this means for our objective function $f$. Having recalled these results we
can immediately clarify the power of such an oracle.

###  2.1 Notions from Complex analysis

Given a sufficiently smooth function
$f:\mathcal{D}\subseteq\mathbb{R}^{n}\to\mathbb{R}$, the most important
results from this section clarify what can be said about the complex extension
of $f$. That is, what kind of structure is preserved if we evaluate $f$ over
some subset of $\mathbb{C}^{n}$ instead of merely over $\mathcal{D}$. To be
clear, when we consider $f:\mathbb{C}^{n}\to\mathbb{C}$, $f:z\mapsto f(z)$,
one usually identifies $\mathbb{C}^{n}$ with $\mathbb{R}^{2n}$ and writes
$f:(x,y)\mapsto f(x+iy)$ with $f(x,y)=u(x,y)+iv(x,y)$ for some functions $u$
and $v$. In what follows we mainly rely on [Kra00, KP02]

A function $f:\mathcal{D}\subseteq\mathbb{R}^{n}\to\mathbb{R}$ is
(real)-analytic on $U\subseteq\mathcal{D}$, denoted $f\in C^{\omega}(U)$, if
for any $x^{\prime}\in U$ there exists a neighbourhood $U^{\prime}\ni
x^{\prime}$ such that for all $x\in U^{\prime}$ the map $f$ can be locally
represented by the convergent power series
$f(x)=\sum_{n=0}^{\infty}a_{n}(x-x^{\prime})^{n}$. This means, for example,
that we have a convergent Taylor series.

Let $\alpha\in\mathbb{N}^{n}$ be a multi-index, that is
$\alpha=(\alpha_{1},\alpha_{2},\dots,\alpha_{n})$,
$|\alpha|=\sum^{n}_{i=1}\alpha_{i}$, $\alpha!=\prod^{n}_{i=1}\alpha_{i}!$ and
most importantly
$\partial_{x}^{\alpha}=\prod_{i=1}^{n}\partial_{x_{i}}^{\alpha_{i}}$. This
allows for a clean description of multivariate Taylor’s series.

######  Theorem 2.1 (Multivariate Taylor’s series [KP02, Section 2.2]).

Let $f\in C^{\omega}(\mathcal{D})$, $\mathcal{D}\subseteq\mathbb{R}^{n}$, then
for each $x_{0}\in\mathcal{D}$ there exists a neighbourhood
$U\subseteq\mathcal{D}$ containing $x_{0}$ such that for each $(x_{0}+h)\in U$

$f(x_{0}+h)=\sum_{|\alpha|\in\mathbb{N}}\frac{1}{\alpha!}(\partial_{x})^{\alpha}f(x_{0})h^{\alpha}.$
(2.1)

To indicate the difference between the real and complex setting, we denote by
$f\in C^{\omega}(\mathcal{D})$ the fact that the real-valued function $f$ is
analytic over the domain $\mathcal{D}\subseteq\mathbb{R}^{n}$, whereas $g\in
H(\Omega)$ denotes that the complex-valued function $g$ is holomorphic
(complex differentiable) over $\Omega\subseteq\mathbb{C}^{n}$.

Although at our core we interested in problems over $\mathbb{R}^{n}$, in
general, however, we will not work over $\mathbb{R}^{n}$, but rather over a
subset of $\mathbb{C}^{n}$. To that end, the following Lemma asserts the
possibility of extending a real-valued domain to a complex-valued domain.

######  Lemma 2.2 (From real-analytic to complex-analytic).

If $f:\mathbb{R}^{n}\to\mathbb{R}$ is real-analytic, then there exists some
open set $\Omega\subseteq\mathbb{C}^{n}$ such that
$\mathbb{R}^{n}\subseteq\Omega$ and $f$ is complex-analytic (or equivalently,
holomorphic) on $\Omega$.

###### Proof.

To be remotely self-contained we recall this basic proof. Since $f\in
C^{\omega}(\mathbb{R}^{n})$ for any $x^{\circ}\in\mathbb{R}^{n}$, there is a
$r>0$ and sequence $(a_{n})_{n\in\mathbb{N}}$ such that the local power series
representation of $f$, that is,
$f(x)=\sum_{n=0}^{\infty}a_{n}(x-x^{\circ})^{n}$ converges for all $x$ in the
open ball $B_{r}(x^{\circ})$. Recall that the radius of convergence $r$ equals
the distance (in $\mathbb{C}^{n}$) from $z^{\circ}:=(x^{\circ},0\cdot
i)\in\mathbb{C}^{n}$ to a singularity, that is, the most nearby point to
$z^{\circ}$ where $f$ fails to be holomorphic, e.g., see [Kra00, Section 2.3].
Hence, the local complex extension follows naturally by considering the same
power series representation via $(a_{n})_{n\in\mathbb{N}}$ yet, for all
$z\in\mathbb{C}^{n}$ in the open ball $B_{r}(z^{\circ})$. The result follows
after applying the same construction for any $x^{\circ}$. ∎

Although perhaps expected, $\Omega$ cannot always be assumed to be of the form
$\Omega\simeq\mathbb{R}^{n}\times(-\epsilon,\epsilon)^{n}$, $\epsilon>0$222For
example, consider $f(x)=\sum^{\infty}_{k=1}a_{k}\frac{1}{1+kx^{2}}$ for some
appropriate sequence $(a_{k})_{k}$. The set of poles contains the set of
points $(1/\sqrt{k})i$, $k\in\mathbb{N}_{>0}$, such that $\Omega$ can never
contain a the full strip of the form
$\mathbb{R}^{n}\times(-\epsilon,\epsilon)^{n}$, $\epsilon>0$.. The set
$\Omega$ is however difficult to quantify. To circumvent this we will mainly
assume that $f:\mathcal{D}\subseteq\mathbb{R}^{n}\to\mathbb{R}$ is holomorphic
over $\mathcal{D}\times(-\bar{\delta},\bar{\delta})$ for some
$\bar{\delta}>0$. This holds for example when $f\in
C^{\omega}(\mathbb{R}^{n})$ and $\mathcal{D}\subset\mathbb{R}^{n}$ is bounded.
One could circumvent this implicit description of $\Omega$ by simply assuming
$f$ to be real-entire333The function $f$ is real-entire when it has a globally
convergent power series representation., but this is a rather strong
assumption, for instance, consider the well-known example function
$f(x)=1/(1+x^{2})$, which fails to be holomorphic at $\pm i$.

The main technical tool is one of the most celebrated tools from Complex
analysis

######  Theorem 2.3 (Multivariate Cauchy-Riemann equations).

Let $f:\Omega\subseteq\mathbb{C}^{n}\to\mathbb{C}$ and define $z:=x+iy$
together with the functions $u$ and $v$ such that the map $f$ is defined by
$x+iy\mapsto f(x+iy)=u(x,y)+iv(x,y)$. Then, if $f$ is holomorphic over
$\Omega$, the Cauchy-Riemann equations

$\partial_{x_{i}}u=\partial_{y_{i}}v,\quad-\partial_{x_{i}}v=\partial_{y_{i}}u,\quad
i=1,2,\dots,n$ (2.2)

are satisfied.

For proofs and more on Theorem 2.3 see [Kra00, KW17]. The converse to Theorem
2.3 is not immediate and usually requires further assumptions [GM78].

We end with an important connection between real-analyticity and strict-
convexity.

######  Lemma 2.4 (Generalized principle of isolated zeroes).

Let $f\in C^{\omega}(\mathbb{R}^{n})$. If $f$ is constant over some open
subset $A\subseteq\mathbb{R}^{n}$ with $\mathrm{vol}(A)>0$, then $f$ is
constant over the whole of $\mathbb{R}^{n}$.

###### Proof.

The proof is a simple consequence of a far reaching result by Łojasiewicz,
e.g., see [Theorem 6.3.3][KP02]. ∎

######  Corollary 2.5 (Real-analytic implies local strict convexity).

For any non-constant convex function $f\in C^{\omega}(\mathbb{R}^{n})$, there
is a line segment $\ell\subset\mathbb{R}^{n}$ such that $f|_{\ell}$ is
strictly convex.

Corollary 2.5 is important to take into account when comparing algorithms, we
do not allow for the most generic convex objective functions.

###  2.2 The imaginary trick

The most straight-forward method of numerical differentiation is the finite-
difference method [Vui+07, Chapter 3] which is based on the approximation

$\partial_{x}f(x)=\frac{f(x+h)-f(x)}{h}+O(h).$

However, as explained before, for sufficiently small $h$ this method suffers
from numerical drawbacks (subtractive cancellation errors). To overcome this
we can use an elegant idea which builds upon some ideas from Complex analysis.
This approach was arguably introduced in [LM67] with the first concrete
complex-step approach appearing in [ST98] and with later elaborations in
[MSA03, ASM15, Abr+18].

Instead of $f$ mapping $\mathbb{R}\to\mathbb{R}$, enlarge its domain to
$\mathbb{C}$, that is, let with some abuse of notation $f(z)=u(z)+iv(z)$ with
$z=x+iy\in\mathbb{C}$ such that $f(x)=f(z)|_{y=0,v(x)=0}$. Now if $f$ is
holomorphic, that is $f\in H(\Omega)$ for $\mathbb{R}\subset\Omega$, we can
appeal to the Cauchy-Riemann equations, telling us that
$\partial_{x}u=\partial_{y}v$ and hence

$\partial_{x}u=\lim_{h\downarrow 0}\frac{v\big{(}x+i(y+h)\big{)}-v(x+iy)}{h}.$
(2.3)

Then since $f(x)=f(z)|_{y=0,v(x)=0}$ and $\Im(f)=v$, we obtain from (2.3)

$\partial_{x}f(x)=\lim_{h\downarrow 0}\frac{\Im\big{(}f(x+ih)\big{)}}{h},$
(2.4)

or non-asymptotically, $\partial_{x}f\approx\Im(f(x+ih))/h$. Hence, there are
no finite-differences, merely a single function evaluation. Moreover, taking
the imaginary part of the Taylor series of $f(x+ih)$ around $x$ (which we can
do due to the holomorphic assumption), we find that

$f(x+ih)=f(x)+\partial_{x}f(x)ih-\frac{1}{2}\partial_{x}^{2}f(x)h^{2}-\frac{1}{6}\partial_{x}^{3}f(x)ih^{3}+O(h^{4})$
(2.5)

and thus

$\partial_{x}f(x)=\frac{\Im\big{(}f(x+ih)\big{)}}{h}+O(h^{2}),$

which is a remarkably powerful observation. Indeed, one can also see the
imaginary trick itself directly from (2.5). Moreover,
$f(x)=\Re(f(x+ih))+O(h^{2})$. So although we do not know $f$, by lifting the
problem to the Complex domain, we can get good approximations of both $f$ and
its gradient using essentially a single function evaluation. Observe that as
for $f(x)=x^{2}$ one has for any $h>0$ that
$\Im(f(x+ih))/h=\partial_{x}f(x)=2x$ such that for locally-quadratic functions
the complex-step approximation is numerically robust.

One can continue this line of thinking and introduce (multi-point) finite
differences again and obtain even smaller errors of the order $O(h^{4})$
[ASM15]. Of course, at the cost of also introducing numerical troubles, so
unless time is expensive, the current complex-step remains favourable.
Moreover, the theory has been extended to the matrix case in [AMH10] and even
to higher-order derivatives [LRD12], plus see [MSA03] for more on the relation
to Automatic Differentiation.

Due to the lack of cancellation errors, the precision of the complex-step
approximation can be effectively arbitrarily high. Especially the aerospace
community has adopted this method actively, the complex-step approximation has
been for example used in air-foil design [GWX17], more impressively it is
reported in [CH04, Page 44] that a value of $h=10^{-100}$ is successfully used
in National Physical Laboratory software. More recently and more related to
our goals, in [NS18] the authors consider using the complex-step coordinate-
wise under complex noise and propose its use in line-search methods. However,
the authors do not provide a complete analysis. The high accuracy potential of
the imaginary trick (2.2) is also mentioned in [SW18, BBN19], but only by
looking at (2.2). We will use the complex-step method in place of the oracle
(1.5) and provide the full analysis.

Especially in the context of optimization algorithms it is desirable to be
able to assert with a selected amount of accuracy that the gradient of the
objective function at the current iterate is sufficiently small, i.e., to
assert some form of local optimality. The next example shows that most
standard numerical approaches to approximate the gradient, but the imaginary
step, fail to posses this desirable property.

(a) $x=-1$

(b) $x=0$

(c) $x=10$

Figure 2.1: Comparison of the different oracles presented in (2.6) on the test
setting of Example 2.6 around a variety of points $x$.

######  Example 2.6 (Numerical oracle quality).

To show the power of the complex-step and concurrently the numerical
difficulties within finite-difference methods we approximate the derivative of
$f(x)=x^{3}$ at $x\in\\{-1,0,10\\}$. We use a forward-difference
($\mathsf{fd}$), central-difference ($\mathsf{cd}$) and complex-step
($\mathsf{cs}$) method

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

and compare the error for $\delta\downarrow 0$. In Figure 2.1 we show the
three different errors as a function of $\delta$. We see that
$f_{\mathsf{cd}}$ and $f_{\mathsf{cs}}$ are of comparable quality, up to some
point, indeed both have an $O(\delta^{2})$ error. However, only the complex-
step can reach machine-precision, whereas the other methods deteriorate
(subtractive cancellation errors become dominant). At this point it is
instrumental to mention that the majority of references use oracles of the
form (2.6a)-(2.6b), hence it is clear that there is room for some numerical
improvements and this indicates the potential of the imaginary trick in the
context of zeroth order optimization.

###  2.3 Lipschitz inequalities

In this part we gather a variety of inequalities which come in useful later.
Note that none of them assumes the underlying function $f$ to be convex.

Using the notation from [Nes03] a function $f$ is said to be an element of
$C^{k,p}_{L}(\mathcal{D})$ when $f$ is $k$ times continuously differentiable
with additionally having its $p^{\mathrm{th}}$-derivative being $L$-Lipschitz
over $\mathcal{D}\subseteq\mathbb{R}^{n}$. That is, if $f\in
C^{1,1}_{L_{1}(f)}(\mathbb{R}^{n})$, then $f$ has a Lipschitz gradient,

$\|\nabla f(x)-\nabla f(y)\|_{2}\leq L_{1}(f)\|x-y\|_{2},\quad\forall
x,y\in\mathbb{R}^{n}$ (2.7)

and indeed, (2.7) implies that

$f(y)\geq f(x)+\langle\nabla f(x),y-x\rangle+\frac{1}{2L_{1}(f)}\|\nabla
f(x)-\nabla f(y)\|_{2}^{2},\quad\forall x,y\in\mathbb{R}^{n}$ (2.8)

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}$,
e.g., see [Zho18]. Moreover, as [Nes11, Equation (6)]

$|f(y)-f(x)-\langle\nabla
f(x),y-x\rangle|\leq\frac{1}{2}L_{1}(f)\|x-y\|_{2}^{2}.$ (2.9)

It follows from [Nes03, Lemma 1.2.4] that if $f\in C^{2}(\mathbb{R}^{n})$ and
satisfies

$\|\nabla^{2}f(x)-\nabla^{2}f(y)\|_{2}\leq L_{2}(f)\|x-y\|_{2}\quad\forall
x,y\in\mathbb{R}^{n},$ (2.10)

for some constant $L_{2}(f)\geq 0$, denoted $f\in
C^{2,2}_{L_{2}}(\mathbb{R}^{n})$, then,

$|f(y)-f(x)-\langle\nabla
f(x),y-x\rangle-\frac{1}{2}\langle\nabla^{2}f(x)(y-x),y-x\rangle|\leq\frac{1}{6}L_{2}(f)\|x-y\|_{2}^{3},\quad\forall
x,y\in\mathbb{R}^{n}.$ (2.11)

See that (2.10) is equivalent to

$|u^{\mathsf{T}}\nabla^{2}f(x)u-u^{\mathsf{T}}\nabla^{2}f(y)u|\leq
L_{2}(f)\|x-y\|_{2}\quad\forall x,y\in\mathbb{R}^{n},\;u\in\mathbb{S}^{n-1},$
(2.12)

which is commonly referred to as $f$ being $3$rd order smooth (cf. [BP16]).
Now it follows directly from the definition of a derivative that $f\in
C^{2,2}_{L_{2}(f)}(\mathbb{R}^{n})$ implies that

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

##  3 An imaginary oracle

In this section we present the main tool of this work. Using ideas from
[Nes11, NS17] we provide the complex-step generalization of the result in
[NY83, FKM04].

To state the main result of this section we define a uniform complex-step
$\delta$-smoothed version of $f$ by

$f_{\delta}(x):=\mathbb{E}_{v\sim\mathbb{B}^{n}}\left[\Re\left(f(x+i\delta
v)\right)\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. It is true
that if one has prior structural knowledge of $f$, one could consider a
different solid to smooth over. Note that if the domain
$\mathcal{D}\subseteq\mathbb{R}^{n}$ of $f$ is for example bounded, then, by
Lemma 2.2, there exists $\bar{\delta}$ such that $f_{\delta}$ is well-defined
for any $\delta\in[0,\bar{\delta})$. To aid readability we will mention it
throughout, but from now on we will assume that Assumption 3.1 holds.

######  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(\bar{\delta},\bar{\delta})^{n}\subset\mathbb{C}^{n}$ for
some $\bar{\delta}\in(0,1)$.

We remark that the $1$ in the interval of $\bar{\delta}$ is larger than
practically needed and merely helps in asserting that for any
$\delta\in(0,\bar{\delta})$ one has $\delta^{p+1}\leq\delta^{p}$ for any
$p\in\mathbb{N}_{\geq 0}$. Also, note that $\mathcal{D}$ need not be bounded,
for example, consider $f(x)=\|Ax-b\|_{2}^{2}$ over $\mathbb{R}^{n}$, which is
a real-entire function and as such $f\in H(\mathbb{C}^{n})$.

######  Lemma 3.2 (Approximation quality of the complex-step function).

Let $f\in C^{\omega}(\mathcal{D})\cap C^{1,1}_{L_{1}(f)}(\mathcal{D})$ admit a
holomorphic extension to
$\mathcal{D}\times(-\bar{\delta},\bar{\delta})^{n}\subset\mathbb{C}^{n}$.
Then, for $f_{\delta}$ as in (3.1) and any $x\in\mathcal{D}$ there exists some
constant $C_{0}>0$ such that for any $\delta\in[0,\bar{\delta})$

$\displaystyle\left|f_{\delta}(x)-f(x)\right|$
$\displaystyle\leq\frac{C_{0}\delta^{2}L_{1}(f)}{2},$ (3.2a)
$\displaystyle\left|f_{\delta}(x)-f(x)\right|$
$\displaystyle\lesssim\frac{\delta^{2}L_{1}(f)}{2}.$ (3.2b)

The bound (3.2a) will be used in the generic case with
$\delta\in[0,\bar{\delta})$, whereas (3.2b) provides for asymptotic (in
$\delta$) insights without universal constants.

###### Proof.

It follows from the definition (3.1),
$\mathbb{E}_{v\sim\mathbb{B}^{n}}[vv^{\mathsf{T}}]=\frac{1}{n+2}\cdot
I_{n}$444To see this, the covariance matrix
$\mathbb{E}_{v\sim\mathbb{B}^{n}}[vv^{\mathsf{T}}]$ must be isotropic. Hence,
for any $i\in[n]$
$\mathbb{E}_{v\sim\mathbb{B}^{n}}v_{i}^{2}=\frac{1}{n}\mathbb{E}_{v\sim\mathbb{B}^{n}}\|v\|_{2}^{2}$
equals $(C/n)\int^{1}_{0}r^{n-1}r^{2}dr=(C/n)\frac{1}{n+2}$, for $r$ the norm
(radius) of $v\in\mathbb{B}^{n}$. To find $C$, use the fact that we integrate
over a probability measure, that is, $1=C\int^{1}_{0}r^{n-1}dr$, which implies
that $(C/n)=1$ from where the result follows. and the Taylor expansion (2.5)
that for all $x\in\mathcal{D}$

$\displaystyle f_{\delta}(x)$
$\displaystyle=\mathbb{E}_{v\sim\mathbb{B}^{n}}\left[f(x)-\frac{1}{2}(\delta
v)^{\mathsf{T}}\nabla^{2}f(x)(\delta v)+O(\|\delta v\|_{2}^{4})\right]$
$\displaystyle=f(x)-\frac{1}{2}\delta^{2}\mathbb{E}_{v\sim\mathbb{B}^{n}}\left[v^{\mathsf{T}}\nabla^{2}f(x)v\right]+O(\delta^{4})$
$\displaystyle=f(x)-\frac{1}{2}\delta^{2}\frac{1}{n+2}\mathrm{Tr}(\nabla^{2}f(x))+O(\delta^{4}).$

Since $\nabla^{2}f(x)\preceq L_{1}(f)\cdot I_{n}$ and $\bar{\delta}<1$, there
exists some constant $C_{0}>0$ such that (3.2a) holds. Similarly, for
sufficiently small $\delta$ the second-order term will dominate which yields
(3.2b). ∎

In contrast to other works (cf. [Nes11]), it is important to remark that
$f_{\delta}$ does not always belongs to the same function class as $f$. For
example, $f(x):=1/(1+x^{2})$ has a Lipschitz continuous gradient with
$L_{1}(f)=2$, while $\lim_{\delta\uparrow 1}f(0+i\delta)=+\infty$ and as such
this makes the function $f_{\delta}$ not necessarily a member of
$C^{1,1}_{L_{1}(f_{\delta})}$. Also, convexity of $f$ does not necessarily
carry over, see Example 3.3 below.

(a) $\delta=10^{-1}$

(b) $\delta=10^{0}$

(c) $\delta=10^{1}$

Figure 3.1: Comparison of the effect of $\delta$ on the non-convexity of
$\Re(f(x+i\delta))$ for $f(x)=x^{4}$.

######  Example 3.3 (Convexity under complex lifting).

Let $f:\mathbb{R}\to\mathbb{R}$ be real-entire, then we find that

$\Re(f(x+i\delta
v))=f(x)+\sum^{\infty}_{k=1}(-1)^{k}\frac{f^{(2k)}(x)}{(2k)!}(\delta v)^{2k}.$

If $f(x)=x^{2}$, then $\Re(f(x+i\delta v))=x^{2}-(\delta v)^{2}$ which is a
convex function in $x$, regardless of the choice of $\delta>0$ and
$v\in\mathbb{R}$. However, for $f(x)=x^{4}$ one finds that $\Re(f(x+i\delta
v))=x^{4}-6x^{2}(\delta v)^{2}+(\delta v)^{4}$, which is a non-convex function
in $x$ for any $\delta>0$ and non-zero $v$, e.g., see that
$\partial_{x}^{2}\Re(f(x+i\delta v))=12x^{2}-12(\delta v)^{2}$ or see Figure
3.1.

As $f_{\delta}$ need not be convex and might not have a Lipschitz continuous
gradient, analyses of algorithms based on the oracle (3.5) require some extra
machinery beyond [NS17].

Now, we state one of the key contributions, which is the integral
representation of $\nabla f_{\delta}$. Lemma 3.4 below is the complex-step
version of the approach as proposed in [NY83, Section 9.3], see e.g., [FKM04,
Lemma 1].

######  Lemma 3.4 (The gradient of the complex-step function).

Let $f\in C^{\omega}(\mathcal{D})$ admit a holomorphic extension to
$\mathcal{D}\times(-\bar{\delta},\bar{\delta})^{n}\subset\mathbb{C}^{n}$, 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})$

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

###### Proof.

Under the assumption that $f$ is holomorphic over a subset of
$\mathbb{C}^{n}$, one can appeal to the Cauchy-Riemann equations [Rud87,
Theorem 11.2] from where it follows that for $n=1$

$\nabla_{x}\int_{-\delta}^{\delta}\Re\left(f(x+iv)\right)dv=\left.\Im\left(f(x+iu)\right)\right|^{\delta}_{-\delta}.$

The additional scaling factor $n/\delta$ is further explained below. For
$n\geq 2$, recall that Stokes’ theorem tells us that
$\int_{\Omega}d\omega=\int_{\partial\Omega}\iota^{\star}\omega$ for
$\iota:\partial\Omega\hookrightarrow\Omega$ [Theorem 7.2.8][AMR88]. After a
combination of Stokes’ theorem and again the Cauchy-Riemann equations we
obtain a result similar to the one as stated in [NY83, FKM04]:

$\nabla_{x}\int_{\delta\mathbb{B}^{n}}\Re\left(f(x+iv)\right)dv=\int_{\delta\mathbb{S}^{n-1}}\Im\left(f(x+iu)\right)\frac{u}{\|u\|_{2}}du.$
(3.4)

Regarding the scaling factor $u/\|u\|_{2}$ see for example [Lee13, Proposition
15.21]. Moreover, due to the assumption that the vectors $v$ and $u$ are
uniformly distributed one obtains

$\displaystyle f_{\delta}(x)=$
$\displaystyle\mathbb{E}_{v\sim\mathbb{B}^{n}}\left[\Re\left(f(x+i\delta
v)\right)\right]=\frac{1}{\mathrm{vol}(\delta\mathbb{B}^{n})}\int_{\delta\mathbb{B}^{n}}\Re\left(f(x+iv)\right)dv,$

and similarly,

$\displaystyle\mathbb{E}_{u\sim\mathbb{S}^{n-1}}\left[\Im\left(f(x+i\delta
u)\right)u\right]=\frac{1}{\mathrm{vol}(\delta\mathbb{S}^{n-1})}\int_{\delta\mathbb{S}^{n-1}}\Im\left(f(x+iu)\right)\frac{u}{\|u\|_{2}}du.$

Combining the aforementioned results with the fact that
$\mathrm{vol}(\delta\mathbb{B}^{n})/\mathrm{vol}(\delta\mathbb{S}^{n-1})=\delta/n$
allows for concluding the proof. ∎

The proof of Lemma 3.4 clearly shows why we went through the Cauchy-Riemann
equations, we critically rely on them to pass from (3.1) to (3.3). Moreover,
Lemma 3.4 provides us with a natural single-point oracle candidate

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

Then, the following result allows for showing consistency below and implies in
particular that
$\mathbb{E}_{u\sim\mathbb{S}^{n-1}}[uu^{\mathsf{T}}]=(1/n)I_{n}$.

######  Lemma 3.5 (Integration over the $n$-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 udu=x.$ (3.6)

###### Proof.

First, rewrite (3.6) as $n\cdot\int_{\mathbb{S}^{n-1}}uu^{\mathsf{T}}du\,x$,
where $uu^{\mathsf{T}}$ represents an outer-product, that is,
$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})$, 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, $f\in C^{r}(\mathcal{D})$ for all $r\geq 1$ and as
such 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$. Now, since

$\langle\nabla f(x),u\rangle=\lim_{\delta\downarrow
0}\frac{1}{\delta}\Im\left(f(x+i\delta u)\right),$

we observe from (3.3) that the approximation is asymptotically consistent,
that is,

$\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 udu\overset{\eqref{equ:int:sphere}}{=}\nabla f(x).$ (3.7)

Exactly this construction was the key observation in [ADX10, Nes11] to reduce
oracle variance and indeed. This observation does not hold for other single-
point oracles (cf. [FKM04, Section 1.1]).

######  Proposition 3.6 (Gradient approximation quality).

Let $f\in C^{\omega}(\mathcal{D})\cap C^{2,2}_{L_{2}(f)}(\mathcal{D})$,
$L_{2}(f)>0$, admit a holomorphic extension to
$\mathcal{D}\times(-\bar{\delta},\bar{\delta})^{n}\subseteq\mathbb{C}^{n}$.
Then, for all $x\in\mathcal{D}$ there is a constant $C_{1}\geq 0$ such that
for any $\delta\in[0,\bar{\delta})$ one has

$\displaystyle\|\nabla f_{\delta}(x)-\nabla f(x)\|_{2}$
$\displaystyle\leq\frac{C_{1}n\delta^{2}L_{2}(f)}{6},$ (3.8a)
$\displaystyle\|\nabla f_{\delta}(x)-\nabla f(x)\|_{2}$
$\displaystyle\lesssim\frac{n\delta^{2}L_{2}(f)}{6}.$ (3.8b)

###### Proof.

Using a similar approach as in [Nes11, Lemma 3] we see that from the Taylor
series expansion of $f(x+i\delta u)$ around $x$ that

$\displaystyle\begin{array}[]{lcl}\|\nabla f_{\delta}(x)-\nabla
f(x)\|_{2}&\overset{\eqref{equ:unif:consistent}}{\leq}&\frac{n}{\delta\cdot\mathrm{vol}(\mathbb{S}^{n-1})}\int_{\mathbb{S}^{n-1}}\left|\Im\left(f(x+i\delta
u)\right)-\delta\langle\nabla f(x),u\rangle\right|\|u\|_{2}du\\\
&\overset{\eqref{equ:third:deriv:bound}}{\leq}&\frac{n}{\delta\cdot\mathrm{vol}(\mathbb{S}^{n-1})}\int_{\mathbb{S}^{n-1}}\left(\frac{\delta^{3}\cdot
L_{2}(f)}{6}+O(\delta^{5})\right)\|u\|_{2}du\\\ &=&\frac{1}{6}\delta^{2}\cdot
n\cdot L_{2}(f)+O(n\cdot\delta^{4}).\end{array}$

Since $\bar{\delta}<1$, this implies that there is a constant $C_{1}\geq 0$
such that (3.8a) holds. For sufficiently small $\delta$ the second order term
$\delta^{2}$ will dominate which yields (3.8b). ∎

We see that our simple single-point approach allows for an error of the form
$O(\delta^{2})$ which is what [Nes11] could get using a Gaussian smoothed
multi-point approach under the assumption that $f\in C^{2,2}_{L_{2}(f)}$.
Removing the constant in (3.8a) Proposition 3.6 is possible by appealing to
for example the Cauchy integral representation theorem, but at the cost of
demanding non-traditional bounds on $f$ over a complex domain.

##  4 Convex optimization

The goal in this section is to solve convex optimization problems inspired by
the zeroth order framework as proposed in [Nes11, NS17], yet, by using the
complex-step oracle $g_{\delta}$ as given by (3.5). The algorithm for the
unconstrained case is detailed in Algorithm 4.1.

1: Input: initial iterate $x_{0}\in\mathbb{R}^{n}$, stepsizes
$\\{\mu_{k}\\}_{k\geq 0}$, smoothing parameters $\\{\delta_{k}\\}_{k\geq 0}$.

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

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

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

5: set $x_{k+1}=x_{k}-\mu_{k}\cdot g_{\delta_{k}}(x_{k})$

6: end for

Algorithm 4.1 Complex $\delta$-smoothed zeroth order unconstrained
optimization

What is more, given some non-empty convex compact set
$\mathcal{K}\subset\mathcal{D}\subseteq\mathbb{R}^{n}$, we also consider
constrained convex optimization problems of the form

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

again, using the zeroth order framework. This algorithm is detailed in
Algorithm 4.2. Here,
$\Pi_{\mathcal{K}}:\mathcal{D}\subseteq\mathbb{R}^{n}\to\mathcal{K}$ denotes
the projection operator.

1: Input: initial iterate $x_{0}\in\mathbb{R}^{n}$, stepsizes
$\\{\mu_{k}\\}_{k\geq 0}$, smoothing parameters $\\{\delta_{k}\\}_{k\geq 0}$.

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

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

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

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

6: end for

Algorithm 4.2 Complex $\delta$-smoothed zeroth order constrained optimization

To characterize the effectiveness of both Algorithm 4.1 and Algorithm 4.2, we
need to bound the variance of the oracle (3.5). We observe the same attractive
property as in [Nes11, NS17], there is no need to assume boundedness of the
second moment of our random oracle.

######  Lemma 4.3 (Oracle variance).

Let $f\in C^{\omega}(\mathcal{D})\cap C^{2,2}_{L_{2}(f)}(\mathcal{D})$ admit a
holomorphic extension to
$\mathcal{D}\times(-\bar{\delta},\bar{\delta})^{n}\subseteq\mathbb{C}^{n}$.
Then, for $g_{\delta}(x)$ as in (3.5) and any $\delta\in[0,\bar{\delta})$
there is a constant $C_{2}\geq 0$ such that

$\displaystyle\mathbb{E}_{u\sim\mathbb{S}^{n-1}}\left[\|g_{\delta}(x)\|_{2}^{2}\right]$
$\displaystyle\leq\frac{C_{2}n^{2}L_{2}(f)^{2}\delta^{4}}{36}+n\|\nabla
f(x)\|_{2}^{2},$ (4.2a)
$\displaystyle\mathbb{E}_{u\sim\mathbb{S}^{n-1}}\left[\|g_{\delta}(x)\|_{2}^{2}\right]$
$\displaystyle\lesssim\frac{n^{2}L_{2}(f)^{2}\delta^{4}}{36}+n\|\nabla
f(x)\|_{2}^{2}.$ (4.2b)

###### Proof.

First, observe from Algorithm 4.1 that

$\mathbb{E}_{u\sim\mathbb{S}^{n-1}}\left[\|g_{\delta}(x)\|_{2}^{2}\right]=\frac{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].$

Then, since $f\in C^{2,2}_{L_{2}}(\mathcal{D})$ one has

$\displaystyle\begin{array}[]{lcl}\Im\left(f(x+i\delta
u)\right)&=&\Im\left(f(x+i\delta u)\right)-\langle\nabla f(x),\delta
u\rangle+\langle\nabla f(x),\delta u\rangle\\\
&\overset{\eqref{equ:third:deriv:bound}}{\leq}&\frac{1}{6}\delta^{3}\cdot
L_{2}(f)+O(\delta^{5})+\langle\nabla f(x),\delta u\rangle,\end{array}$

such that

$\mathbb{E}_{u\sim\mathbb{S}^{n-1}}\left[\|g_{\delta}(x)\|_{2}^{2}\right]\overset{\eqref{equ:int:sphere}}{\leq}\frac{1}{36}n^{2}L_{2}(f)^{2}\delta^{4}+n\|\nabla
f(x)\|_{2}^{2}+O(n^{2}\cdot\delta^{6}).$

Hence, since $\bar{\delta}<1$ there is constant $C_{2}\geq 0$ such that for
any $\delta\in[0,\bar{\delta})$ (4.2a) holds. For sufficiently small $\delta$,
the lowest order term in $\delta$ will dominate and hence we recover (4.2b). ∎

Lemma 4.3 follows from the fact that the leading term in (2.5) is $O(h^{3})$.
This should be contrasted with standard (purely real) results where this term
would be $O(h^{2})$. Note, [Nes11, Theorem 4] obtains a bound similar to
(4.2a) under the assumption that $f\in C^{2,2}_{L_{2}(f)}$, yet, using a
multi-point scheme. Also note that a variety of other zeroth order schemes
have $O(1/\delta^{p})$ showing up in their variance bounds which is however
due to the measurement noise, i.e., they have only access to $f(x)+\xi$, which
is not the setting considered in this paper.

######  Theorem 4.4 (Convergence rate of Algorithm 4.1 and Algorithm 4.2).

Let $f:\mathcal{D}\subseteq\mathbb{R}^{n}\to\mathbb{R}$ be a convex function
that admits a holomorphic extension to
$\mathcal{D}\times(-\bar{\delta},\bar{\delta})^{n}\subseteq\mathbb{C}^{n}$ for
some $0<\bar{\delta}<1$ and let $\mathcal{K}\subset\mathcal{D}$ be a compact
convex set. Suppose that $f$ has a Lipschitz gradient and Hessian, that is,
(2.7) and (2.10) hold, for non-zero constants $L_{1}(f)$ and $L_{2}(f)$,
respectively. Let $\\{x_{k}\\}_{k\geq 0}$ be the sequence of iterates
generated by either Algorithm 4.1 or Algorithm 4.2, with constant stepsize
$\mu_{k}\equiv\mu=1/({2nL_{1}(f)})$ and smoothing parameters
$\\{\delta_{k}\\}_{k\geq 0}\subset(0,\bar{\delta})$. Then, there are constants
$C_{1},C_{2},\dots,C_{b},C_{c}$ such that for any $K\geq 1$,

  1. (i)

If $\delta_{k}\in(0,\bar{\delta})$, then the sequence
$\\{x_{k}\\}_{k=0}^{K-1}$ satisifes

$\displaystyle\frac{1}{K}\sum_{k=0}^{K-1}\mathbb{E}\left[f(x_{k})-f(x^{\star})\right]\leq$
$\displaystyle\frac{2nL_{1}(f)}{K}\|x_{0}-x^{\star}\|_{2}^{2}+\frac{C_{1}nL_{2}(f)}{K}\|x_{0}-x^{\star}\|_{2}\sum_{k=0}^{K-1}\delta_{k}^{2}$
(4.3) $\displaystyle+\frac{nL_{2}(f)^{2}}{2L_{1}(f)K}$
$\displaystyle\left(C_{2}\left(\sum_{k=0}^{K-1}\delta_{k}^{2}\right)^{2}+C_{3}\sum_{k=0}^{K-1}\delta_{k}^{2}\left(\sum_{k=0}^{K-1}\delta_{k}^{4}\right)^{\frac{1}{2}}+C_{4}\sum_{k=0}^{K-1}\delta_{k}^{4}\right);$

  2. (ii)

If $\delta_{k}\equiv\delta$ for some $\delta\in(0,\bar{\delta})$ then
$\bar{x}_{K-1}:=\frac{1}{K}\sum_{k=0}^{K-1}x_{k}$ satisfies

$\mathbb{E}\left[f(\bar{x}_{K-1})-f(x^{\star})\right]\leq\frac{2nL_{1}(f)}{K}\|x_{0}-x^{\star}\|_{2}^{2}+{C_{a}nL_{2}(f)\delta^{2}}\|x_{0}-x^{\star}\|_{2}+\frac{C_{b}nL_{2}(f)^{2}K\delta^{4}}{L_{1}(f)};$
(4.4)

  3. (iii)

While if $\delta_{k}=\min\\{\bar{\delta},\delta/(k+1)\\}$, for some
$\delta\in(0,\bar{\delta})$, then

$\mathbb{E}\left[f(\bar{x}_{K-1})-f(x^{\star})\right]\leq\frac{n}{K}\left(\sqrt{2L_{1}(f)}\|x_{0}-x^{\star}\|_{2}+\delta^{2}\frac{L_{2}(f)C_{c}}{\sqrt{2L_{1}(f)}}\right)^{2}.$
(4.5)

###### Proof.

We do the proof for Algorithm 4.1 and show that the extension to Algorithm 4.2
is immediate. The abbreviations $\mathrm{C.S.}$ and $\mathrm{T.I.}$ denote the
Cauchy-Schwarz and triangle inequality, respectively.

Let $x^{\star}$ be any optimal solution. Define
$r_{k}:=\|x_{k}-x^{\star}\|_{2}$ for $k\geq 0$. Then, there are non-negative
constants $C_{1},C_{2}$ such that

$\displaystyle\begin{array}[]{lcl}\mathbb{E}_{u_{k}\sim\mathbb{S}^{n-1}}\left[r_{k+1}^{2}|x_{k}\right]&=&\mathbb{E}_{u_{k}\sim\mathbb{S}^{n-1}}\left[r_{k}^{2}-2\mu_{k}\langle
g_{\delta_{k}}(x_{k}),x_{k}-x^{\star}\rangle+\mu_{k}^{2}\|g_{\delta_{k}}(x_{k})\|_{2}^{2}|x_{k}\right]\\\
&\overset{\eqref{equ:unif:grad}}{=}&r_{k}^{2}-2\mu_{k}\langle\nabla
f_{\delta_{k}}(x_{k}),x_{k}-x^{\star}\rangle+\mu_{k}^{2}\mathbb{E}_{u_{k}\sim\mathbb{S}^{n-1}}\left[\|g_{\delta_{k}}(x_{k})\|_{2}^{2}|x_{k}\right]\\\
&\overset{\mathrm{C.S.}}{\leq}&r_{k}^{2}-2\mu_{k}\langle\nabla
f(x_{k}),x_{k}-x^{\star}\rangle+2\mu_{k}\|\nabla f_{\delta_{k}}(x_{k})-\nabla
f(x_{k})\|_{2}\|x_{k}-x^{\star}\|_{2}\\\
&&+\mu_{k}^{2}\mathbb{E}_{u_{k}\sim\mathbb{S}^{n-1}}\left[\|g_{\delta_{k}}(x_{k})\|_{2}^{2}|x_{k}\right]\\\
&\overset{\eqref{equ:grad:delta:approx:error},\eqref{equ:g:var}}{\leq}&r_{k}^{2}-2\mu_{k}\left(f(x_{k})-f(x^{\star})\right)+\mu_{k}\frac{C_{1}n\delta_{k}^{2}L_{2}(f)}{3}\|x_{k}-x^{\star}\|_{2}\\\
&&+\mu_{k}^{2}\left(\frac{1}{36}C_{2}n^{2}L_{2}(f)^{2}\delta_{k}^{4}+n\|\nabla
f(x_{k})\|_{2}^{2}\right)\\\
&\overset{\eqref{equ:grad:Lipschitz:bound2}}{\leq}&r_{k}^{2}-2\mu_{k}\left(f(x_{k})-f(x^{\star})\right)+\mu_{k}\frac{C_{1}n\delta_{k}^{2}L_{2}(f)}{3}\|x_{k}-x^{\star}\|_{2}\\\
&&+\mu_{k}^{2}\left(\frac{1}{36}C_{2}n^{2}L_{2}(f)^{2}\delta_{k}^{4}+n2L_{1}(f)(f(x_{k})-f(x^{\star}))\right),\end{array}$

which, upon rearranging and taking expectation on $u_{1},\dots,u_{k-1}$,
yields

$\displaystyle\mathbb{E}\left[f(x_{k})-f(x^{\star})\right]$
$\displaystyle\leq\frac{1}{\mu}\left(\mathbb{E}\left[r_{k}^{2}\right]-\mathbb{E}\left[r_{k+1}^{2}\right]\right)+\frac{\delta_{k}^{2}C_{1}nL_{2}(f)}{3}\mathbb{E}\left[r_{k}\right]+\frac{\mu
C_{2}n^{2}L_{2}(f)^{2}\delta_{k}^{4}}{36}$
$\displaystyle\leq\frac{1}{\mu}\left(\mathbb{E}\left[r_{k}^{2}\right]-\mathbb{E}\left[r_{k+1}^{2}\right]\right)+\frac{\delta_{k}^{2}C_{1}nL_{2}(f)}{3}\sqrt{\mathbb{E}\left[r_{k}^{2}\right]}+\frac{\mu
C_{2}n^{2}L_{2}(f)^{2}\delta_{k}^{4}}{36}.$

Then, summing the above inequalities for $k$ from $0$ to $K-1$ gives

$\displaystyle\sum_{k=0}^{K-1}\mathbb{E}\left[f(x_{k})-f(x^{\star})\right]\leq\frac{1}{\mu}\left(r_{0}^{2}-\mathbb{E}\left[r_{K}^{2}\right]\right)+\frac{C_{1}nL_{2}(f)}{3}\sum_{k=0}^{K-1}\delta_{k}^{2}\sqrt{\mathbb{E}\left[r_{k}^{2}\right]}+\frac{\mu
C_{2}n^{2}L_{2}(f)^{2}}{36}\sum_{k=0}^{K-1}\delta_{k}^{4}.$ (4.6)

For $k\geq 0$, letting

$s_{k}=\sqrt{\mathbb{E}\left[r_{k}^{2}\right]},$

we then have

$s_{K}^{2}\overset{\eqref{ineq:f_k_suboptimality_sum}}{\leq}s_{0}^{2}+\frac{\mu
C_{1}nL_{2}(f)}{3}\sum_{k=0}^{K-1}\delta_{k}^{2}s_{k}+\frac{\mu^{2}C_{2}n^{2}L_{2}(f)^{2}}{36}\sum_{k=0}^{K-1}\delta_{k}^{4}.$
(4.7)

Applying Lemma A.1 to (4.7) yields

$\displaystyle s_{K}$ $\displaystyle\leq\frac{\mu
C_{1}nL_{2}(f)}{6}\sum_{k=0}^{K-1}\delta_{k}^{2}+\left(s_{0}^{2}+\frac{\mu^{2}C_{2}n^{2}L_{2}(f)^{2}}{36}\sum_{k=0}^{K-1}\delta_{k}^{4}+\left(\frac{\mu
C_{1}nL_{2}(f)}{6}\sum_{k=0}^{K-1}\delta_{k}^{2}\right)^{2}\right)^{\frac{1}{2}}$
$\displaystyle\overset{\mathrm{T.I.}}{\leq}\frac{\mu
C_{1}nL_{2}(f)}{3}\sum_{k=0}^{K-1}\delta_{k}^{2}+s_{0}+\left(\frac{\mu^{2}C_{2}n^{2}L_{2}(f)^{2}}{36}\sum_{k=0}^{K-1}\delta_{k}^{4}\right)^{\frac{1}{2}}.$
(4.8)

Substituting (4.8) into (4.6), we get (4.3) from

$\begin{split}&\,\frac{1}{K}\sum_{k=0}^{K-1}\mathbb{E}\left[f(x_{k})-f(x^{\star})\right]\\\
\leq&\,\frac{1}{\mu
K}\|x_{0}-x^{\star}\|_{2}^{2}+\frac{C_{1}nL_{2}(f)}{3K}\sum_{k=0}^{K-1}\delta_{k}^{2}\left(\frac{\mu
C_{1}nL_{2}(f)}{3}\sum_{k=0}^{K-1}\delta_{k}^{2}+s_{0}+\left(\frac{\mu^{2}C_{2}n^{2}L_{2}(f)^{2}}{36}\sum_{k=0}^{K-1}\delta_{k}^{4}\right)^{\frac{1}{2}}\right)\\\
&+\frac{\mu C_{2}n^{2}L_{2}(f)^{2}}{36K}\sum_{k=0}^{K-1}\delta_{k}^{4}\\\
\leq&\,\frac{1}{\mu
K}\|x_{0}-x^{\star}\|_{2}^{2}+\frac{C_{1}nL_{2}(f)}{3K}\|x_{0}-x^{\star}\|_{2}\sum_{k=0}^{K-1}\delta_{k}^{2}+\frac{\mu
C_{1}^{2}n^{2}L_{2}(f)^{2}}{9K}\left(\sum_{k=0}^{K-1}\delta_{k}^{2}\right)^{2}\\\
&\,\qquad+\frac{\mu
C_{1}C_{2}^{\frac{1}{2}}n^{2}L_{2}(f)^{2}}{18K}\left(\sum_{k=0}^{K-1}\delta_{k}^{2}\right)\left(\sum_{k=0}^{K-1}\delta_{k}^{4}\right)^{\frac{1}{2}}+\frac{\mu
C_{2}n^{2}L_{2}(f)^{2}}{36K}\sum_{k=0}^{K-1}\delta_{k}^{4}.\end{split}$ (4.9)

Without loss of generality, let both $C_{1}\geq 1$ and $C_{2}\geq 1$. Then, if
$\delta_{k}\equiv\delta\in(0,\bar{\delta})$, we derive (4.4) from inequality
(4.9)

$\displaystyle\,\mathbb{E}\left[f(\bar{x}_{K-1})-f(x^{\star})\right]\leq\frac{1}{K}\sum_{k=0}^{K-1}\mathbb{E}\left[f(x_{k})-f(x^{\star})\right]$
$\displaystyle\leq$ $\displaystyle\,\frac{1}{\mu
K}\|x_{0}-x^{\star}\|_{2}^{2}+\frac{C_{1}nL_{2}(f)\delta^{2}}{3}\|x_{0}-x^{\star}\|_{2}+\frac{\mu
C_{2}n^{2}L_{2}(f)^{2}\delta^{4}}{36}$ $\displaystyle+\frac{\mu
C_{1}^{2}n^{2}L_{2}(f)^{2}K\delta^{4}}{9}+\frac{\mu
C_{1}C_{2}^{\frac{1}{2}}n^{2}L_{2}(f)^{2}\sqrt{K}\delta^{4}}{18}$
$\displaystyle\leq$ $\displaystyle\,\frac{1}{\mu
K}\|x_{0}-x^{\star}\|_{2}^{2}+\frac{C_{1}nL_{2}(f)\delta^{2}}{3}\|x_{0}-x^{\star}\|_{2}+\frac{\mu
C_{1}^{2}C_{2}n^{2}L_{2}(f)^{2}K\delta^{4}}{4}$ $\displaystyle=$
$\displaystyle\,\frac{2nL_{1}(f)}{K}\|x_{0}-x^{\star}\|_{2}^{2}+\frac{C_{1}nL_{2}(f)\delta^{2}}{3}\|x_{0}-x^{\star}\|_{2}+\frac{C_{1}^{2}C_{2}nL_{2}(f)^{2}K\delta^{4}}{8L_{1}(f)}.$

If however $\delta_{k}=\min\left\\{\bar{\delta},\delta/(k+1)\right\\}$ for all
$k\geq 0$, we derive (4.5) from inequality (4.9), the zeta function
inequalities (A.1) and the convexity of $f$ via

$\displaystyle\mathbb{E}\left[f(\bar{x}_{K-1})-f(x^{\star})\right]\leq$
$\displaystyle\frac{1}{K}\sum_{k=0}^{K-1}\mathbb{E}\left[f(x_{k})-f(x^{\star})\right]$
$\displaystyle\leq$ $\displaystyle\frac{1}{\mu
K}\|x_{0}-x^{\star}\|_{2}^{2}+\delta^{2}\frac{\pi^{2}C_{1}nL_{2}(f)}{18K}\|x_{0}-x^{\star}\|_{2}+\delta^{4}\frac{\mu\pi^{4}C_{1}^{2}n^{2}L_{2}(f)^{2}}{324K}$
$\displaystyle+\delta^{4}\frac{\mu\pi^{4}C_{1}C_{2}^{\frac{1}{2}}n^{2}L_{2}(f)^{2}}{108\sqrt{90}K}+\delta^{4}\frac{\mu\pi^{4}C_{2}n^{2}L_{2}(f)^{2}}{3240K}$
$\displaystyle=$ $\displaystyle\,\frac{1}{\mu
K}\|x_{0}-x^{\star}\|_{2}^{2}+\delta^{2}\frac{nL_{2}(f)C_{3}}{K}\|x_{0}-x^{\star}\|_{2}+\delta^{4}\frac{\mu
n^{2}L_{2}(f)^{2}C_{4}}{K}$ $\displaystyle=$
$\displaystyle\,\frac{2nL_{1}(f)}{K}\|x_{0}-x^{\star}\|_{2}^{2}+\delta^{2}\frac{nL_{2}(f)C_{3}}{K}\|x_{0}-x^{\star}\|_{2}+\delta^{4}\frac{nL_{2}(f)^{2}C_{4}}{2L_{1}(f)K}$
$\displaystyle\leq$
$\displaystyle\,\frac{n}{K}\left(\sqrt{2L_{1}(f)}\|x_{0}-x^{\star}\|_{2}+\delta^{2}\frac{L_{2}(f)C_{5}}{\sqrt{2L_{1}(f)}}\right)^{2}.$

This completes the proof for Algorithm 4.1.

Now consider Algorithm 4.2. Let $r_{k}=\|x_{k}-x^{\star}\|_{2}$, then since
$\mathcal{K}$ is a compact convex set it follows that

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

Hence, we can appeal to proof of Algorithm 4.1 and conclude on the convergence
rate. ∎

Combining all asymptotic (in $\delta$) error bounds (3.2b), (3.8b), (4.2b)
with (4.5), then all universal constants can be taken to be $1$ and as such it
can be observed that if $\delta$ is sufficiently small indeed, then, one needs

$O\left(\frac{nL_{1}(f)\|x_{0}-x^{\star}\|_{2}^{2}}{\epsilon}\right)$ (4.10)

iterations to have
$\mathbb{E}\left[f(\bar{x}_{K-1})-f(x^{\star})\right]\leq\epsilon$. This is
precisely the first-order, i.e., with $\nabla f(x)$ being available,
complexity scaled by the dimension $n$ [Nes03, Section 2.1.5].

##  5 Strongly convex optimization

In this part the previous results on unconstrained optimization are extended
to the setting of $f$ 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+\frac{\tau(f)}{2}\|y-x\|_{2}^{2},\quad\forall
x,y\in\mathcal{D}.$ (5.1)

In particular (5.1) implies

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

If additionally $f\in C^{1,1}_{L_{1}(f)}$, then

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

######  Theorem 5.1 (Convergence rate of Algorithm 4.1 under strong
convexity).

Let $f:\mathcal{D}\subseteq\mathbb{R}^{n}\to\mathbb{R}$ be a
$\tau(f)$-strongly convex function that admits a holomorphic extension to
$\mathcal{D}\times(-\bar{\delta},\bar{\delta})^{n}\subseteq\mathbb{C}^{n}$ for
some $0<\bar{\delta}<1$. Suppose that $f$ has a Lipschitz gradient and
Hessian, that is, (2.7) and (2.10) hold, for non-zero constants $L_{1}(f)$ and
$L_{2}(f)$, respectively. Let $\\{x_{k}\\}_{k\geq 0}$ be the sequence of
iterates generated by Algorithm 4.1 with constant stepsize
$\mu_{k}\equiv\mu={1}/({2nL_{1}(f)})$ and the sequence of smoothing parameters
defined for all $k\geq 0$ by $\delta_{k}=\min\\{\bar{\delta},\delta/(k+1)\\}$
for some $\delta\in(0,\bar{\delta})$. Then, there are non-negative constants
$C_{a},C_{b}$ such that if we let

$\nu_{\tau,\delta}=\frac{\delta^{2}nL_{1}(f)}{\tau(f)}\left(\|x_{0}-x^{\star}\|_{2}^{2}\frac{C_{a}L_{2}(f)}{L_{1}(f)}+\delta^{2}\frac{C_{b}L_{2}(f)^{2}}{L_{1}(f)^{2}}\right).$
(5.4)

then, for any $K\geq 0$, $x_{K+1}$ satisfies

$\mathbb{E}[f(x_{K+1})-f(x^{\star})]\leq\frac{1}{2}L_{1}(f)\left(\nu_{\tau,\delta}+\left(1-\frac{\tau(f)}{4nL_{1}(f)}\right)^{K}(\|x_{0}-x^{\star}\|_{2}^{2}-\nu_{\tau,\delta})\right).$
(5.5)

###### Proof.

Let $x^{\star}$ be any optimal solution. Define
$r_{k}:=\|x_{k}-x^{\star}\|_{2}$ and
$\rho_{k+1}:=\mathbb{E}_{u_{k}\sim\mathbb{S}^{n-1}}\left[r_{k+1}^{2}|x_{k}\right]$
for $k\geq 0$. Then, it follows from the proof of Theorem 4.4 that there are
non-negative constants $C_{1},C_{2}$ such that

$\displaystyle\begin{array}[]{lcl}\rho_{k+1}&\overset{\eqref{equ:grad:Lipschitz:bound2}}{\leq}&r_{k}^{2}-2\mu_{k}\left(f(x_{k})-f(x^{\star})\right)+\mu_{k}\frac{C_{1}n\delta_{k}^{2}L_{2}(f)}{3}\|x_{k}-x^{\star}\|_{2}\\\
&&+\mu_{k}^{2}\left(\frac{1}{36}C_{2}n^{2}L_{2}(f)^{2}\delta_{k}^{4}+n2L_{1}(f)(f(x_{k})-f(x^{\star}))\right)\\\
&\overset{\eqref{equ:strong:ineq}}{\leq}&\left(1-\mu_{k}\frac{\tau(f)}{2}\right)r_{k}^{2}+\mu_{k}\frac{C_{1}n\delta_{k}^{2}L_{2}(f)}{3}r_{k}+\mu_{k}^{2}\left(\frac{1}{36}C_{2}n^{2}L_{2}(f)^{2}\delta_{k}^{4}\right).\end{array}$

Therefore

$\displaystyle\rho_{k+1}$
$\displaystyle\leq\left(1-\frac{\tau(f)}{4nL_{1}(f)}\right)\rho_{k}+\delta_{k}^{2}\frac{C_{1}L_{2}(f)}{6L_{1}(f)}\sqrt{\rho_{k}}+\delta_{k}^{4}\frac{C_{2}L_{2}(f)^{2}}{144L_{1}(f)^{2}}$
(5.6a)
$\displaystyle\overset{\eqref{equ:Lip:strongly:convex}}{\leq}\rho_{k}+\delta_{k}^{2}\frac{C_{1}L_{2}(f)}{6L_{1}(f)}\sqrt{\rho_{k}}+\delta_{k}^{4}\frac{C_{2}L_{2}(f)^{2}}{144L_{1}(f)^{2}}$
(5.6b)

Now summing up over $k$ and using a telescoping argument as in the proof of
Theorem 4.4 yields

$\rho_{K}\leq\rho_{0}+\sum^{K-1}_{k=0}\delta_{k}^{2}\frac{C_{1}L_{2}(f)}{6L_{1}(f)}\sqrt{\rho_{k}}+\sum^{K-1}_{k=0}\delta_{k}^{4}\frac{C_{2}L_{2}(f)^{2}}{144L_{1}(f)^{2}}.$
(5.7)

Then, applying Lemma A.1 to (5.7) yields

$\displaystyle\sqrt{\rho_{K}}$
$\displaystyle\leq\frac{1}{2}\sum^{K-1}_{k=0}\delta_{k}^{2}\frac{C_{1}L_{2}(f)}{6L_{1}(f)}+\left(\rho_{0}^{2}+\sum^{K-1}_{k=0}\delta_{k}^{4}\frac{C_{2}L_{2}(f)^{2}}{144L_{1}(f)^{2}}+\left(\frac{1}{2}\sum^{K-1}_{k=0}\delta_{k}^{2}\frac{C_{1}L_{2}(f)}{6L_{1}(f)}\right)^{2}\right)^{\frac{1}{2}}$
$\displaystyle\overset{\mathrm{T.I.}}{\leq}\sum^{K-1}_{k=0}\delta_{k}^{2}\frac{C_{1}L_{2}(f)}{6L_{1}(f)}+\rho_{0}+\left(\sum^{K-1}_{k=0}\delta_{k}^{4}\frac{C_{2}L_{2}(f)^{2}}{144L_{1}(f)^{2}}\right)^{\frac{1}{2}}.$

Plugging this bound on $\sqrt{\rho_{K}}$ back into (5.6a) yields

$\displaystyle\rho_{K+1}\leq$
$\displaystyle\left(1-\frac{\tau(f)}{4nL_{1}(f)}\right)\rho_{K}+\delta_{K}^{4}\frac{C_{2}L_{2}(f)^{2}}{144L_{1}(f)^{2}}$
$\displaystyle+\delta_{K}^{2}\frac{C_{1}L_{2}(f)}{6L_{1}(f)}\left(\sum^{K-1}_{k=0}\delta_{k}^{2}\frac{C_{1}L_{2}(f)}{6L_{1}(f)}+\rho_{0}+\left(\sum^{K-1}_{k=0}\delta_{k}^{4}\frac{C_{2}L_{2}(f)^{2}}{144L_{1}(f)^{2}}\right)^{\frac{1}{2}}\right).$

Recall the zeta bound (A.1) and let us select
$\delta_{k}=\min\\{\bar{\delta},\delta/(k+1)\\}$ such that

$\displaystyle\rho_{K+1}\leq$
$\displaystyle\left(1-\frac{\tau(f)}{4nL_{1}(f)}\right)\rho_{K}+\delta^{4}\frac{\pi^{2}C_{1}^{2}L_{2}(f)^{2}}{216L_{1}(f)^{2}}+\delta^{2}\rho_{0}\frac{C_{1}L_{2}(f)}{6L_{1}(f)}+\delta^{4}\frac{C_{2}L_{2}(f)^{2}}{144L_{1}(f)^{2}}+\delta^{4}\frac{\pi^{2}C_{1}C_{2}^{\frac{1}{2}}L_{2}(f)^{2}}{72\sqrt{90}L_{1}(f)^{2}}.$

Hence, there are constants $C_{3}$ and $C_{4}$ such that

$\rho_{K+1}\leq\left(1-\frac{\tau(f)}{4nL_{1}(f)}\right)\rho_{K}+\delta^{2}\rho_{0}C_{3}\frac{L_{2}(f)}{L_{1}(f)}+\delta^{4}C_{4}\frac{L_{2}(f)^{2}}{L_{1}(f)^{2}}.$

Then, for $\nu_{\tau,\delta}$ as in (5.4) we get

$(\rho_{K+1}-\nu_{\tau,\delta})\leq\left(1-\frac{\tau(f)}{4nL_{1}(f)}\right)(\rho_{K}-\nu_{\tau,\delta})\leq\left(1-\frac{\tau(f)}{4nL_{1}(f)}\right)^{K}(\rho_{0}-\nu_{\tau,\delta}).$

At last, it follows from (2.9) that
$\mathbb{E}[f(x_{k})-f(x^{\star})]\leq\frac{1}{2}L_{1}(f)\rho_{k}$, which
concludes the proof. ∎

Now, contrasting Theorem 4.4 with Theorem 5.1 it follows directly from (5.4)
and (5.5) that to have
$\mathbb{E}\left[f({x}_{K+1})-f(x^{\star})\right]\leq\epsilon$ one needs a
sufficiently small $\delta$ satisfying

$\delta\leq O\left(\sqrt{\frac{\epsilon\tau(f)}{nL_{1}(f)^{2}}}\right)$ (5.8)

and

$O\left(\frac{nL_{1}(f)}{\tau(f)}\log\left(\frac{L_{1}(f)\|x_{0}-x^{\star}\|_{2}^{2}}{\epsilon}\right)\right)$
(5.9)

iterations. Comparing the rate (5.9) to (4.10), then we see — as is known, see
e.g., [Nes11]— that strong convexity allows for a faster scheme.

##  6 Non-convex optimization

At last we consider a local minima in a possibly non-convex program. Our line
of proof is different from [Nes11, Section 7] as our smoothed function
$f_{\delta}$ does not necessarily has a Lipshitz continuous gradient. In this
particular setting, the set $\mathcal{D}$ functions as an open neighbourhood
of a local minima.

######  Theorem 6.1 (Convergence rate of Algorithm 4.1 to a local minima).

Let $f:\mathcal{D}\subseteq\mathbb{R}^{n}\to\mathbb{R}$ be a — not necessarily
convex — function that admits a holomorphic extension to
$\mathcal{D}\times(-\bar{\delta},\bar{\delta})^{n}\subseteq\mathbb{C}^{n}$ for
some $0<\bar{\delta}<1$. Suppose that $f$ has a Lipschitz gradient and
Hessian, that is, (2.7) and (2.10) hold, for non-zero constants $L_{1}(f)$ and
$L_{2}(f)$, respectively. Let $\\{x_{k}\\}_{k\geq 0}$ be the sequence of
iterates generated by Algorithm 4.1 with constant stepsize
$\mu_{k}\equiv\mu={1}/({nL_{1}(f)})$ and the sequence of smoothing parameters
defined for all $k\geq 0$ by $\delta_{k}=\min\\{\bar{\delta},\delta/(k+1)\\}$
for some $\delta\in(0,\bar{\delta})$. Let $x^{\star}\in\mathcal{D}$ be a local
minima of $f$, then, there is a non-negative constant $C_{a}$ such that for
all $K\geq 1$ the sequence $\\{x_{k}\\}_{k=0}^{K-1}$ satisfies

$\displaystyle\frac{1}{K}\sum^{K-1}_{k=0}\mathbb{E}\left[\|\nabla
f(x_{k})\|_{2}^{2}\right]\leq$
$\displaystyle\frac{n}{K}\Big{[}2L_{1}(f)\left(f(x_{0})-f(x^{\star})\right)$
(6.1)
$\displaystyle+\delta^{2}L_{2}(f)\left(\delta^{2}L_{2}(f)(n+1)C_{a}+\sqrt{2nL_{1}(f)(f(x_{0})-f(x^{\star}))}\right)\Big{]}.$

###### Proof.

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+\frac{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+\frac{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.8a) and (4.2a) 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}\frac{C_{1}n\delta_{k}^{2}L_{2}(f)}{6}\|\nabla
f(x_{k})\|_{2}$ $\displaystyle+\frac{1}{2}\mu_{k}^{2}L_{1}(f)\left(n\|\nabla
f(x_{k})\|_{2}^{2}+\frac{C_{2}n^{2}L_{2}(f)^{2}\delta_{k}^{4}}{36}\right).$

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

$\displaystyle\frac{\mu}{2}\mathbb{E}\left[\|\nabla f(x_{k})\|_{2}\right]^{2}$
$\displaystyle\leq\frac{\mu}{2}\mathbb{E}\left[\|\nabla
f(x_{k})\|_{2}^{2}\right]$ (6.2)
$\displaystyle\leq\mathbb{E}[f(x_{k})-f(x_{k+1})]+\frac{C_{1}\delta_{k}^{2}L_{2}(f)}{6L_{1}(f)}\mathbb{E}\left[\|\nabla
f(x_{k})\|_{2}\right]+\frac{C_{2}L_{2}(f)^{2}\delta_{k}^{4}}{72L_{1}(f)}.$

As we consider a local minima, assume that $f(x)\geq f(x^{\star})$ and define
$\phi_{k}:=\mathbb{E}[\|\nabla f(x_{k})\|_{2}]$, then, a telescoping argument
yields

$\displaystyle\phi_{K}^{2}\leq
2nL_{1}(f)(f(x_{0})-f(x^{\star}))+\frac{nC_{1}L_{2}(f)}{3}\sum^{K}_{k=0}\delta_{k}^{2}\phi_{k}+\frac{nC_{2}L_{2}(f)^{2}}{36}\sum^{K}_{k=0}\delta_{k}^{4}.$

Since $\|\nabla f(x_{0})\|_{2}^{2}\leq 2L_{1}(f)|f(x_{0})-f(x^{\star})|$ we
can appeal to Lemma A.1 and obtain

$\displaystyle\phi_{K}\leq\frac{nC_{1}L_{2}(f)}{6}\sum^{K}_{k=0}\delta_{k}^{2}+\left(2nL_{1}(f(x_{0})-f(x^{\star}))+\frac{nC_{2}L_{2}(f)^{2}}{36}\sum^{K}_{k=0}\delta^{4}_{k}+\left(\frac{nC_{1}L_{2}(f)}{6}\sum^{K}_{k=0}\delta_{k}^{2}\right)^{2}\right)^{\frac{1}{2}}.$

Now we consider $\delta_{k}=\min\\{\bar{\delta},\delta/(k+1)\\}$ such that
after using (A.1) and an application of the triangle inequality we obtain

$\displaystyle\phi_{K}$
$\displaystyle\leq\sqrt{2nL_{1}(f)(f(x_{0})-f(x^{\star}))}+\delta^{2}\frac{nC_{1}L_{2}(f)\pi^{2}}{18}+\delta^{2}\frac{\sqrt{n}C_{2}^{\frac{1}{2}}L_{2}(f)\pi^{2}}{6\sqrt{90}}$
$\displaystyle\leq\sqrt{2nL_{1}(f)(f(x_{0})-f(x^{\star}))}+\delta^{2}\frac{nC_{3}L_{2}(f)\pi^{2}}{18}.$

Then, define $\theta_{k}^{2}:=\mathbb{E}\left[\|\nabla
f(x_{k})\|_{2}^{2}\right]$, plug $\phi_{K}$ back into (6.2) and obtain after a
telescoping argument

$\displaystyle\frac{1}{K+1}\sum^{K}_{k=0}\theta^{2}_{k}\leq$
$\displaystyle\frac{n}{K+1}\Big{[}2L_{1}(f)\left(f(x_{0})-f(x^{\star})\right)+\delta^{4}C_{4}L_{2}(f)^{2}$
$\displaystyle+\delta^{2}L_{2}(f)\left(\delta^{2}nL_{2}(f)C_{5}+\sqrt{2nL_{1}(f)(f(x_{0})-f(x^{\star}))}\right)\Big{]}.$

∎

Given (6.1) we see that if $\delta$ is chosen sufficiently small, then, to
achieve $(1/K)\sum^{K-1}_{k=0}\mathbb{E}[\|\nabla
f(x_{k})\|_{2}^{2}]\leq\epsilon$ we need

$K\geq
O\left(\frac{nL_{1}(f)\left(f(x_{0})-f(x^{\star})\right)}{\epsilon}\right).$

##  7 Numerical experiments

Especially in the zeroth order setting as considered here, one has usually no
direct knowledge of the sharpest Lipschitz constants related to the objective
function $f$. In such a case one needs to resort to an approximation and in
fact for the Algorithms to work, the estimated Lipschitz constant should not
under-approximate the real Lipschitz constant. However, as remarked throughout
the literature — for deterministic oracles (cf. [Nes11, Equations (58,59)])— a
larger Lipschitz constant implies a smaller smoothing parameter $\delta$ to
obtain a certain suboptimality gap $\epsilon$, see also (5.8). Hence, having
the freedom to pick $\delta$ arbitrarily small is preferable.

In the following experiments we will show that multi-point oracles do not
allow for this option whereas our single-point method does. At the same time,
if $\delta$ is sufficiently large such that multi-point methods work, we show
that our method achieves a similar or faster convergence rate. All results
above show that the (asymptotic) rates are the same as the ones provided in
[NS17], yet we observe a slightly better dependency on $n$, that is, $n$
versus $2(n+4)$. A full investigation of why our method leads to empirically
faster convergence is left to future work.

As mentioned throughout the introduction, the oracle (1.4) has severe
limitations — most notably very slow convergence — which is why we omitted a
comparison against it, this would not be informative. Instead we compare
against the Gaussian smoothing (GS) oracles from [NS17, Equation (30)], namely

$\displaystyle g_{\delta,\mathsf{fd}}(x)$ $\displaystyle=\frac{f(x+\delta
u)-f(x)}{\delta}u,\quad u\sim\mathcal{N}(0,I_{n}),$ (7.1a) $\displaystyle
g_{\delta,\mathsf{cd}}(x)$ $\displaystyle=\frac{f(x+\delta u)-f(x-\delta
u)}{2\delta}u,\quad u\sim\mathcal{N}(0,I_{n}).$ (7.1b)

As indicated, these oracles are based on finite-difference methods from
numerical differentiation (2.6). When simulating (7.1a) and (7.1b) we will
always use the stepsize $\mu_{k}\equiv 1/(4(n+4)L_{1}(f))$ as given by [NS17,
Equation (55)].

ia Suboptimality gap $f(\bar{x}_{K})-f^{\star}$ for the test function (7.2).

iia Cost $f(\bar{x}_{K})$ for the Huber-like function (7.3).

iiia Cost $f(\bar{x}_{K})$ for logistic regression (7.5).

ib Suboptimality gap $f({x}_{K})-f^{\star}$ for the test function (7.2).

iib Cost $f({x}_{K})$ for the Huber-like function (7.3).

iiib Cost $f({x}_{K})$ for logistic regression (7.5).

Figure 7.1: The single-point Complex-smoothing ($\mathrm{CS}$) method
(Algorithm 4.1) compared to the multi-point Gaussian smoothing ($\mathrm{GS}$)
method from [NS17, Equation (54), oracle (7.1a)] on a variety of objective
functions for a time-invariant smoothing parameter $\delta\equiv\delta_{k}$.

###  7.1 Unconstrained convex optimization

Throughout this section, we always compare against oracle (7.1a) unless
indicated differently.

######  Example 7.1 (Quadratic test function).

First, we compare Algorithm 4.1 to [Nes11, NS17] on an ill-conditioned version
of the test function555Better known as the “worst function in the world”. from
[Nes03, Section 2.1.2]

$f_{n}(x)={L}\left(\frac{1}{2}\left[(x^{(1)})^{2}+\sum^{n-1}_{i=1}(x^{(i+1)}-x^{(i)})^{2}+(x^{(n)})^{2}\right]-x^{(1)}\right)$
(7.2)

for $x_{0}=0$, $L=10^{-8}$, $L_{1}(f)=4L$ and $(x^{\star})^{(i)}=1-i/(n+1)$
with $x^{(i)}$ denoting the $i^{\mathrm{th}}$ component of
$x\in\mathbb{R}^{n}$. Indeed, on such a quadratic function, the complex-step
and two-point oracle as proposed in [Nes11, NS17] have a theoretically
equivalent performance. However, fixing $n=5$ and letting the smoothing
parameter approach machine precision shows the critical difference between a
one-point and multi-point method. See Figure 1ia and Figure 1ib for a single
realization, where we compare the performance of $\bar{x}_{K}$ and $x_{K}$,
respectively. Especially the performance of $x_{K}$ is significantly better
under the complex-step algorithm.

######  Example 7.2 (Smooth approximate $\ell_{1}$-regularization).

Following [FG16], we are interested in solving a smoothly approximated version
of a $\ell_{1}$-regularized convex program. Specifically, we consider the
pseudo-Huber loss given by

$\psi_{\mu}(x)=\mu\sum^{m}_{i=1}\left(\sqrt{1+\frac{x_{i}^{2}}{\mu^{2}}}-1\right)$
(7.3)

and are interested in minimizing the objective

$f(x)=\frac{1}{2}\|Ax-b\|_{2}^{2}+\lambda\psi_{\mu}(x)$ (7.4)

over $x\in\mathbb{R}^{n}$ for some $\lambda>0$ and data
$A\in\mathbb{R}^{m\times n}$, $b\in\mathbb{R}^{m}$. Now, it follows from
[FG16, Lemma 2] that $L_{1}(f)=\lambda/\mu+\|A^{\top}A\|_{2}$, whereas it
follows from [FG16, Lemma 6] that $L_{2}(f)=\lambda/\mu^{2}$. Now, we again
compare Algorithm 4.1 to the method proposed in [Nes11, NS17]. Here we let $A$
and $b$ be random with unit covariance matrices for $m=4$, $n=2$. Moreover,
$\lambda=\mu=10^{-4}$. We show the costs $f(\bar{x}_{K})$ and $f(x_{k})$ in
Figure 1iia and Figure 1iib, respectively, for a decreasing smoothing
parameter $\delta$. Again, a difference in numerical stability can be
observed.

######  Example 7.3 (Logistic regression).

For the function

$f(x)=\frac{1}{m}\sum^{m}_{i=1}\log\left(1+\exp\left(-y_{i}a_{i}^{\mathsf{T}}x\right)\right)$
(7.5)

with $y_{i}\in\\{-1,1\\}$ and $a_{i}\sim\mathcal{N}(0,I_{n})$, one can derive
that $L_{1}(f)=\frac{1}{4}\|A\|_{2}$, for $A$ the matrix filled with the
vectors $a_{i}$. Let $m=100$ and $n=2$, then again, we can observe mild
numerical issues in Figure 1iiia for a sufficiently small $\delta$ if one uses
the averaged iterate $\bar{x}_{K}$. However, these are significantly more
pronounced when one is interested in a faster algorithm and uses $x_{k}$ as
shown in Figure 1iiib.

(a) Suboptimality gap $f(\bar{x}_{K})-f^{\star}$ for
$f(x)=\frac{1}{2}\|x\|_{2}^{2}$ with $n=1$.

(b) Suboptimality gap $f(\bar{x}_{K})-f^{\star}$ for
$f(x)=\frac{1}{2}\|x\|_{2}^{2}$ with $n=100$.

(c) Suboptimality gap $f(\bar{x}_{K})-f^{\star}$ for
$f(x)=\frac{1}{2}\|x\|_{2}^{2}$ with $n=10000$.

(d) Suboptimality gap $f({x}_{K})-f^{\star}$ for
$f(x)=\frac{1}{2}\|x\|_{2}^{2}$ with $n=1$.

(e) Suboptimality gap $f({x}_{K})-f^{\star}$ for
$f(x)=\frac{1}{2}\|x\|_{2}^{2}$ with $n=100$.

(f) Suboptimality gap $f({x}_{K})-f^{\star}$ for
$f(x)=\frac{1}{2}\|x\|_{2}^{2}$ with $n=10000$.

Figure 7.2: The single-point Complex-smoothing ($\mathrm{CS}$) method
(Algorithm 4.1) compared to the multi-point Gaussian smoothing ($\mathrm{GS}$)
method from [NS17, Equation (55)] for oracle (7.1a) (GSfd) and oracle (7.1b)
(GScd) as a function of the dimension $n$ for a time-invariant smoothing
parameter $\delta\equiv\delta_{k}$.

######  Example 7.4 (The effect of dimension and strong convexity).

Figure 7.1 showed the effect of numerical cancellation, but also a clear
difference in convergence speed. Looking back at Example 2.6 and comparing
oracles (7.1a) and (7.1b) to the complex-step oracle (3.5) we expect that a
speed-up should be visible be passing from the forward-difference to the
central-difference oracle (7.1b). In Figure 7.2 we compare these oracles under
an increase in dimension $n$ on the function $f(x)=\frac{1}{2}\|x\|_{2}^{2}$
for $x_{0}=(1/\sqrt{n})\mathbbold{1}_{n}$. We observe that the severe
numerical issues prevail, yet slightly less for higher dimensions.
Concurrently, we observe throughout that averaging is more robust, but
significantly slower. What is more, although the central-difference oracle
(7.1b) does provide for a speed-up compared to (7.1a), the convergence remains
slower than Algorithm 4.1 under the complex-step oracle (3.5). Also, plugging
in $\tau(f)=1$, $L_{1}(f)=1$ and $\|x_{0}\|_{2}^{2}=1$ into the provided
strong convexity rate from Theorem 5.1 we observe that $(1/2)(1-1/(4n))^{K}$
does capture the convergence rate of $f(x_{k})$ to $f(x^{\star})$. We compare
this against the rate $(1/2)(1-1/(8(n+4)))^{K}$ provided by [NS17, Equation
(57)].

###  7.2 Constrained convex optimization

Next we do a control example where we iteratively update the input policy. In
the unconstrained case one could consider the policy iteration schemes
provided in [Faz+18, Mal+19].

######  Example 7.5 (Policy iteration).

We are interested in solving the following Model Predictive Control (MPC)
problem. There, given some prediction horizon $T\geq 1$ one wants to solve

$\displaystyle\operatorname*{\text{minimize}}_{u:=\\{u_{t}\\}_{t=0}^{T-1}\subset\mathbb{R}^{n_{u}}}$
$\displaystyle\sum^{T-1}_{t=0}\langle Qx_{t},x_{t}\rangle+\langle
Ru_{t},u_{t}\rangle+\langle Qx_{T},x_{T}\rangle$ (7.6)
$\displaystyle\operatorname*{\text{subject to}}$ $\displaystyle
x_{t+1}=Ax_{t}+Bu_{t},\quad x_{0}=x^{\prime},$
$\displaystyle\|u_{t}\|_{\infty}\leq 1.$

for the policy
$u^{\star}(x^{\prime})=u_{0}^{\star}(x^{\prime}),\dots,u_{T-1}^{\star}(x^{\prime})$.
Given $u^{\star}(x^{\prime})$, only the first input, that is
$u_{0}^{\star}(x^{\prime})$ is implemented. This input maps $x_{0}$ to
$x_{1}=Ax_{0}+Bu_{0}^{\star}(x_{0})$. Now we set $x^{\prime}=x_{1}$ and
continue this iterative process for as long as we would like the control
horizon to be. In this example we assume that the cost matrices
$Q\in\mathcal{S}^{n_{x}}_{\succeq 0}$ and $R\in\mathcal{S}^{n_{u}}_{\succ 0}$
are given. However, we do not have access to the model parameters
$A\in\mathbb{R}^{{n_{x}}\times n_{x}}$ and $B\in\mathbb{R}^{n_{x}\times
n_{u}}$ directly, we can only simulate the performance. That is, given some
policy $u^{\prime}=\\{u^{\prime}_{t}\\}_{t=0}^{T-1}$ we can simulate the
dynamics and hence evaluate the corresponding cost in (7.6). Here, we let the
tuple $(A,B,Q,R)$ be the standard 2 dimensional Yalmip [Löf04] MPC
instance666Available at https://yalmip.github.io/example/standardmpc/ with
$x_{0}=(3,1)$. We will construct the optimal sequence of control inputs
$u^{\star}_{0}(x_{0}),u^{\star}_{0}(x_{1}),\dots,$ via Yalmip, and denote this
as ${u}^{\star}$. Note that Problem (7.6) is of the form
$\min_{u\in\mathcal{K}}f(u)$ for $\mathcal{K}$ a compact and convex set. Now
we will show the performance of Algorithm 4.2 compared to oracle (7.1b) under
a time-invariant $\delta_{k}\equiv\delta=10^{-16}$ for $K=10^{6}$ and $T=4$.
We use $u_{0}=0\in\mathbb{R}^{n_{u}\cdot T}$ as our initialization and use
$L_{1}(f)=4\cdot 10^{4}$ as our Lipschitz estimate. Note, this is a very rough
estimate based on assumed knowledge of operator norms of $A$ and $B$. We
simulate the MPC scheme for a control horizon of length $15$. In Figure 3ia
and Figure 3ib we show the typical convergence for solving one instance of
(7.6). It should be remarked that the oscillations as seen in Figure 3ib are
due to the optimizer living on the boundary of the feasible set. Then,
iteratively applying the MPC scheme, in Figure 3ic and Figure 3id we show that
both the resulting control policy as well as the controlled trajectory induce
by the complex-step method are close(r) to the optimal one.

###  7.3 Non-convex optimization

We end with a classical example.

######  Example 7.6 (Rosenbrock function).

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

$\operatorname*{\text{minimize}}_{x\in\sqrt{2}\mathbb{B}^{2}}\quad(1-x^{(1)})^{2}+100\left((x^{(2)}-(x^{(1)})^{2}\right)^{2}.$
(7.7)

The minimizing argument of (7.7) is $x^{\star}=(1,1)$, which happens to be the
global optimizer over $\mathbb{R}^{2}$ as well. We will compare Algorithm 4.1
with the stepsize $\mu_{k}\equiv\mu=1/(nL_{1}(f))$ from Theorem 6.1 against
[NS17, Section 7]. Here we use oracle (7.1b) but remark that using oracle
(7.1a) results in effectively the same performance on this problem. We start
both algorithms from the same 4 initial conditions in $\sqrt{2}\mathbb{B}^{2}$
and show the convergence of the suboptimality gap in Figure 3iia and the paths
taken by the algorithms in Figure 3iib. The complex-step method convergences
significantly faster. We could speed-up this method slightly by further
decreasing $\delta=10^{-10}$, however, then the Gaussian smoothing method
would break down.

ia Suboptimality gap $f(\bar{u}_{K})-\widehat{f}^{\star}$ for (7.6).

ic Optimal control policies for Example 7.5.

iia Suboptimality gap $f({x}_{K})-f^{\star}$ for (7.7).

ib Suboptimality gap $f({u}_{K})-\widehat{f}^{\star}$ for (7.6).

id Controlled trajectories resulting from the policies above.

iib Paths taken corresponding to Figure 3iia.

Figure 7.3: The single-point Complex-smoothing ($\mathrm{CS}$) method applied
to (i) a constrained problem (Example 7.5) using Algorithm 4.2 and to (ii) a
non-convex problem (Example 7.6) using Theorem 6.1.

Many more examples are possible, for example, as sketched in [CSV09, Section
1.2] one could consider tuning the regularization parameter in ridge
regression. Comparing to for example the numerical experiment from [NG21] is
less insightful as their smoothing parameter relies on the variance of the
noise.

In conclusion, we see that the performance highlighted in Example 2.6
extrapolates to zeroth order optimization algorithms using those oracles. Here
we remark that if one uses the averaged iterate
$\bar{x}_{K-1}=(1/K)\sum^{K-1}_{k=0}x_{k}$, then, these numerical problems are
less pronounced. This comes of course at the cost of slower convergence. Hence
the complex-step method provides an attractive alternative if a provably fast
and numerically stable algorithm is desired.

##  8 Future work

As highlighted in [AMH10], smoothness is sufficient but not necessary for the
complex-step approximation to work. Moreover, the use of generic dual numbers
— which are closely related to $i=\sqrt{-i}$ — as used in automatic
differentiation might bring about new insights. It would also be interesting
to extend the theory to the class of weakly-convex functions and to
investigate the multi-batch and online settings.

##  Appendix A Appendix

######  Lemma A.1 ([SRB11, Lemma 1]).

Let $\\{t_{k}\\}_{k\geq 0}$ be a non-negative sequence satisfying the
recursion

$t_{K}^{2}\leq T_{K}+\sum_{k=0}^{K-m}\nu_{k}t_{k}\quad\forall K\geq 0,$

where $m\in\\{0,1\\}$ and $\\{T_{K}\\}_{K\geq 0}$ is a non-decreasing sequence
with $T_{0}\geq t_{0}^{2}$ and $\nu_{k}\geq 0$ $\forall k\geq 0$. Then,

$t_{K}\leq\frac{1}{2}\sum_{k=0}^{K-m}\nu_{k}+\left(T_{K}+\left(\frac{1}{2}\sum_{k=0}^{K-m}\nu_{k}\right)^{2}\right)^{\frac{1}{2}}\quad\forall
K\geq 0.$

###### Proof.

Up to a change in summation indices, the proof follows directly from [SRB11,
Lemma 1]. ∎

Throughout a variety of the proofs we appeal to the following inequalities

$\sum^{J}_{j=1}\frac{1}{j^{2}}\leq\frac{\pi^{2}}{6},\quad\sum^{J}_{j=1}\frac{1}{j^{4}}\leq\frac{\pi^{4}}{90},\quad\forall
J\in\mathbb{N},$ (A.1)

which are truncations of the corresponding Riemann zeta functions.

#### Acknowledgements

WJ and DK are supported by the Swiss National Science Foundation under the
NCCR Automation, grant agreement 51NF40_180545. MCY is supported by the Hong
Kong Research Grants Council under the grant 25302420. WJ likes to thank Prof.
Nemirovski for the supplied reference material and Prof. Faulwasser for the
introduction to the imaginary trick.

### 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 _COLT_ , 2010, pp. 28–40
  * [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
  * [AMR88] R Abraham, J.E. Marsden and T Ratiu  “Manifolds, Tensor analysis, and Applications, Second Edition”, Applied Mathematical Sciences  Springer, 1988
  * [APT20] Arya Akhavan, Massimiliano Pontil and Alexandre B Tsybakov  “Exploiting higher order smoothness in derivative-free optimization and continuous bandits”, 2020  arXiv:2006.07862
  * [ASM15] Rafael Abreu, Daniel Stich and Jose Morales  “The Complex-Step-Finite-Difference method”  In _Geophysical Journal International_ 202.1, 2015, pp. 72–93
  * [BBN19] Albert S. Berahas, Richard H. Byrd and Jorge Nocedal  “Derivative-Free Optimization of Noisy Functions via Quasi-Newton Methods”  In _SIAM Journal on Optimization_ 29.2, 2019, pp. 965–993
  * [BCS19] Albert S Berahas, Liyuan Cao and Katya Scheinberg  “Global convergence rate analysis of a generic line search algorithm with noise”, 2019  arXiv:1910.04055
  * [BG18] Krishnakumar Balasubramanian and Saeed Ghadimi  “Zeroth-order nonconvex stochastic optimization: Handling constraints, high-dimensionality and saddle-points”, 2018  arXiv:1809.06474
  * [BM13] Francis Bach and Eric Moulines  “Non-strongly-convex smooth stochastic approximation with convergence rate O(1/n)”  In _Advances in Neural Information Processing Systems_ 26 Curran Associates, Inc., 2013, pp. 773–781
  * [BP16] Francis Bach and Vianney Perchet  “Highly-smooth zero-th order online optimization”  In _Conference on Learning Theory_ , 2016, pp. 1–27
  * [Cai+21] HanQin Cai, Daniel Mckenzie, Wotao Yin and Zhenliang Zhang  “A One-bit, Comparison-Based Gradient Estimator”, 2021  arXiv:2010.02479
  * [CH04] MG Cox and PM Harris  “Software Support for Metrology Best Practice Guide No. 11”, 2004
  * [CSV09] Andrew R. Conn, Katya Scheinberg and Luis N. Vicente  “Introduction to Derivative-Free Optimization”  Society for IndustrialApplied Mathematics, 2009
  * [d’A08] Alexandre d’Aspremont  “Smooth optimization with approximate gradient”  In _SIAM Journal on Optimization_ 19.3 SIAM, 2008, pp. 1171–1183
  * [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
  * [Faz+18] Maryam Fazel, Rong Ge, Sham Kakade and Mehran Mesbahi  “Global Convergence of Policy Gradient Methods for the Linear Quadratic Regulator”  In _Proceedings of the 35th International Conference on Machine Learning_ 80, Proceedings of Machine Learning Research  Stockholmsmässan, Stockholm Sweden: PMLR, 2018, pp. 1467–1476
  * [FG16] Kimon Fountoulakis and Jacek Gondzio  “A second-order method for strongly convex $\ell_{1}$-regularization problems”  In _Mathematical Programming_ 156.1 Springer, 2016, pp. 189–219
  * [FKM04] Abraham Flaxman, Adam Tauman Kalai and H. Brendan McMahan  “Online convex optimization in the bandit setting: gradient descent without a gradient”  In _CoRR_ , 2004
  * [Gas+17] Alexander V Gasnikov et al.  “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
  * [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
  * [GM78] J. D. Gray and S. A. Morris  “When is a Function that Satisfies the Cauchy-Riemann Equations Analytic?”  In _The American Mathematical Monthly_ 85.4 Mathematical Association of America, 1978, pp. 246–256
  * [Gol+19] Daniel Golovin et al.  “Gradientless descent: High-dimensional zeroth-order optimization”, 2019  arXiv:1911.06317
  * [GWX17] Yisheng Gao, Yizhao Wu and Jian Xia  “Automatic differentiation based discrete adjoint method for aerodynamic design optimization on unstructured meshes”  In _Chinese Journal of Aeronautics_ 30.2, 2017, pp. 611–627
  * [HL14] Elad Hazan and Kfir Y Levy  “Bandit Convex Optimization: Towards Tight Bounds.”  In _NIPS_ , 2014, pp. 784–792
  * [Ji+19] Kaiyi Ji, Zhe Wang, Yi Zhou and Yingbin Liang  “Improved zeroth-order variance reduced algorithms and analysis for nonconvex optimization”  In _International Conference on Machine Learning_ , 2019, pp. 3100–3109  PMLR
  * [KP02] Steven G. Krantz and Harold R. Parks  “A Primer of Real Analytic Functions”  Birkhäuser, 2002
  * [Kra00] Steven G. Krantz  “Function Theory of Several Complex Variables”  AMS Chelsea Publishing, 2000
  * [KW17] Jacob Korevaar and Jan Wiegerinck  “Several complex variables”, 2017  URL: https://staff.science.uva.nl/j.j.o.o.wiegerinck/edu/scv/scvboek.pdf
  * [LBM20] Jiaxiang Li, Krishnakumar Balasubramanian and Shiqian Ma  “Stochastic Zeroth-order Riemannian Derivative Estimation and Optimization”, 2020  arXiv:2003.11238
  * [Lee13] John M. Lee  “Introduction to Smooth Manifolds”, Graduate Texts in Mathematics  Springer, 2013
  * [Lia+16] Xiangru Lian et al.  “A Comprehensive Linear Speedup Analysis for Asynchronous Stochastic Parallel Optimization from Zeroth-Order to First-Order”  In _Proceedings of the 30th International Conference on Neural Information Processing Systems_ , 2016, pp. 3062–3070
  * [Liu+20] S. Liu et al.  “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 _arXiv preprint arXiv:1904.11585_ , 2019
  * [Löf04] J. Löfberg  “YALMIP : A Toolbox for Modeling and Optimization in MATLAB”  In _In Proceedings of the CACSD Conference_ , 2004
  * [LRD12] Gregory Lantoine, Ryan P Russell and Thierry Dargent  “Using multicomplex variables for automatic computation of high-order derivatives”  In _ACM Transactions on Mathematical Software_ 38.3, 2012, pp. 1–21
  * [Mal+19] Dhruv Malik et al.  “Derivative-free methods for policy optimization: Guarantees for linear quadratic systems”  In _The 22nd International Conference on Artificial Intelligence and Statistics_ , 2019, pp. 2916–2925  PMLR
  * [MSA03] Joaquim R. R. A. Martins, Peter Sturdza and Juan J. Alonso  “The Complex-Step Derivative Approximation”  In _ACM Trans. Math. Softw._ 29.3 New York, NY, USA: Association for Computing Machinery, 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
  * [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.”, 1983
  * [Ove01] Michael L Overton  “Numerical computing with IEEE floating point arithmetic”  SIAM, 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
  * [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
  * [Rud87] Walter Rudin  “Real and Complex Analysis”, Third edition  McGraw-Hill Education, 1987
  * [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
  * [SMG13] Sebastian U Stich, Christian L Muller and Bernd Gartner  “Optimization of convex functions with random pursuit”  In _SIAM Journal on Optimization_ 23.2 SIAM, 2013, pp. 1284–1309
  * [SRB11] Mark Schmidt, Nicolas Roux and Francis Bach  “Convergence rates of inexact proximal-gradient methods for convex optimization”  In _Advances in neural information processing systems_ 24, 2011, pp. 1458–1466
  * [ST98] William Squire and George Trapp  “Using Complex Variables to Estimate Derivatives of Real Functions”  In _SIAM Review_ 40.1, 1998, pp. 110–112
  * [SV15] Yaron Singer and Jan Vondrák  “Information-theoretic lower bounds for convex optimization with erroneous oracles”  In _Advances in Neural Information Processing Systems_ 28, 2015, pp. 3204–3212
  * [SW18] Jan A Snyman and Daniel N Wilke  “Practical mathematical optimization: basic optimization theory and gradient-based algorithms”  Springer, 2018
  * [Vui+07] Cornelis Vuik, Fredericus Johannes Vermolen, Martin Bastiaan Gijzen and MJ Vuik  “Numerical methods for ordinary differential equations”  VSSD, 2007
  * [Zha+20] Yan Zhang, Yi Zhou, Kaiyi Ji and Michael M Zavlanos  “Improving the Convergence Rate of One-Point Zeroth-Order Optimization using Residual Feedback”, 2020  arXiv:2006.10820
  * [Zho18] Xingyu Zhou  “On the Fenchel Duality between Strong Convexity and Lipschitz Continuous Gradient”, 2018