OE_^ OmOE_^ OmOE_^ Omm!OE_^ mmmmOE_^

# A Gradient Smoothed Functional Algorithm with Truncated Cauchy Random
Perturbations for Stochastic Optimization

Akash Mondal
Indian Institute of Science
<EMAIL_ADDRESS>Prashanth L. A.
Indian Institute of Technology Madras
<EMAIL_ADDRESS>Shalabh Bhatnagar
Indian Institute of Science
<EMAIL_ADDRESS>
###### Abstract

In this paper, we present a stochastic gradient algorithm for minimizing a
smooth objective function that is an expectation over noisy cost samples, and
only the latter are observed for any given parameter. Our algorithm employs a
gradient estimation scheme with random perturbations, which are formed using
the truncated Cauchy distribution from the $\delta$ sphere. We analyze the
bias and variance of the proposed gradient estimator. Our algorithm is found
to be particularly useful in the case when the objective function is non-
convex, and the parameter dimension is high. From an asymptotic convergence
analysis, we establish that our algorithm converges almost surely to the set
of stationary points of the objective function and obtain the asymptotic
convergence rate. We also show that our algorithm avoids unstable equilibria,
implying convergence to local minima. Further, we perform a non-asymptotic
convergence analysis of our algorithm. In particular, we establish here a non-
asymptotic bound for finding an $\epsilon$-stationary point of the non-convex
objective function. Finally, we demonstrate numerically through simulations
that the performance of our algorithm outperforms GSF, SPSA and RDSA by a
significant margin over a few non-convex settings and further validate its
performance over convex (noisy) objectives.

##  1 Introduction

In this paper, we consider the following stochastic optimization (SO) problem:

$\textrm{ Find
}f^{*}:=\inf_{x\in\mathbb{R}^{d}}\left\\{f(x)=\int_{\Xi}F(x,\xi)dP(\xi)\right\\},$
(1)

where $f\colon\mathbb{R}^{d}\to\mathbb{R}$ is a smooth function that could be
highly nonlinear with multiple local minima, $\xi$ is the noise random
variable (r.v.) with support $\Xi$, and $F(x,\xi)$ is a noisy observation of
the function value $f(x)$. We do not assume precise gradient information is
available. Instead, gradients need to be estimated using the aforementioned
noisy observations at certain parameter values.

Stochastic approximation (SA) is an important technique for solving SO
problems. Robbins and Monro (1951) first proposed the SA approach for the
problem of root finding and Kiefer and Wolfowitz (1952) presented the first
application of SA for solving SO problems. Many popular incremental-update
procedures for root finding involving noisy function observations are SA
algorithms. These algorithms are employed for back-propagation in neural
networks (Gawthrop and Sbarbaro, 1990), for solving least squares objectives
(Yao et al., 2009), and finding optimal policies in reinforcement learning
problems (Bertsekas, 2019). As a result, advancements in general SA
methodology have an impact on a wide range of applications.

We consider in this paper a stochastic gradient (SG) algorithm based on
gradient estimates obtained from noisy function observations. SG algorithms
that update without direct (though possibly noisy) observations of the
function gradient are referred to as zeroth order SO algorithms. The first
such algorithm was presented in Kiefer and Wolfowitz (1952), and requires $2d$
function observations per iteration since it perturbs each parameter component
separately.

Such an approach does not scale well in terms of the computational complexity
as the parameter dimension $d$ is increased.

The random directions stochastic approximation (RDSA) method has been
presented in Kushner and Clark (1978b) and also explored in Chin (1997). The
idea here is to randomly perturb all parameter components simultaneously using
random vectors that are uniformly distributed over the surface of the unit
sphere. In Prashanth et al. (2017), independent, symmetric, uniformly
distributed perturbations have been explored and both gradient and Newton RDSA
algorithms have been proposed and analyzed for their asymptotic convergence
properties and rates.

The smoothed functional (SF) algorithm based on independent Gaussian random
perturbations has been independently studied in Katkovnik and Kulchitsky
(1972), Kreimer and Rubinstein (1972), Nesterov and Spokoiny (2017), Bhatnagar
et al. (2013). We shall refer to this algorithm as GSF. The idea underlying
GSF is to approximate the convolution of the objective-function gradient with
a multi-variate Gaussian PDF with the convolution of the objective function
with a scaled Gaussian. Thus, this procedure works with only one simulation
regardless of the parameter dimension. A two-simulation finite-difference
version of the same with lower bias has been studied in Styblinski and Tang
(1990) and Chin (1997). In Bhatnagar (2007), Newton-based algorithms with
biased gradient and Hessian estimates obtained from Gaussian perturbations
have been analyzed for their asymptotic convergence.

The simultaneous perturbation stochastic approximation (SPSA) algorithm, see
Spall (1992), has also gathered attention due to the low computational effort
required in this scheme as well as the ease of implementation. This algorithm
randomly perturbs all parameter components simultaneously by using
perturbation random variates whose properties are commonly satisfied by
independent, symmetric, zero-mean, Bernoulli r.v.

##### Our contributions.

We now summarize our contributions below.
(a) SG with Cauchy perturbations: For solving the SO problem 1, we propose an
SG algorithm, which performs gradient estimation using SF estimates based on
$d$-dimensional truncated Cauchy perturbations inside a sphear of radius
$\delta$.
(b) Asymptotic convergence: We prove that our algorithm converges
asymptotically to the set of local minima of the objective function $f$. Note
here that one can ordinarily prove convergence of an SG algorithm to the set
of equilibria of the associated gradient ODE. These points however also
include local maxima and saddle points (in addition to local minima) of the
given objective function. We however prove by verifying a result from Pemantle
(1990) that the convergence of our algorithm is only to local minima. In fact,
the other fixed points (that are not local minima) are unstable equilibria of
the ODE and which we show are avoided by the algorithm in the limit. A result
of this nature has not been claimed for GSF, to the best of our knowledge.
(c) Asymptotic convergence rate: Our algorithm provides better asymptotic
mean-squared error (AMSE) as compared to SPSA and GSF, which are two very
popular schemes for gradient estimation in the zeroth-order SO context that we
consider. AMSE is a standard metric for comparing the (asymptotic) convergence
rate of different algorthms, cf. Chin (1997), Prashanth et al. (2017), and a
better AMSE is beneficial in simulation optimization applications, where each
function measurement is assumed to be computationally intensive.
(d) Non-asymptotic bound: We also provide non-asymptotic bounds in the spirit
of Ghadimi and Lan (2013), i.e., to an $\epsilon$-solution (see Definition 2
below) of the SO problem mentioned above. In the latter work, the authors
consider a GSF scheme for gradient estimation, and provide an
$O\left(\frac{1}{\epsilon^{2}}\right)$ bound on the number of iterations to
find an $\epsilon$-solution under the assumption that the function $F$
accounting for the noisy observations is smooth. We provide a matching bound
for our proposed algorithm. Further, unlike Ghadimi and Lan (2013), we also
provide non-asymptotic bounds in the case when $F$ is not assumed to be
smooth.
(e) Simulation experiments: Numerical results using a quadratic function,
Rastrigin’s function and Rosenbrock’s function establish that our algorithm
outperforms GSF, SPSA and RDSA algorithm.

##### Comparison to related works.

In Chin (1997), Spall (1992), Bhatnagar and Borkar (2003), Prashanth et al.
(2017), the authors employ random perturbations based gradient estimates
within an SG framework, and mainly show asymptotic convergence to a stationary
point for their algorithms. Our algorithm, on the other hand, is shown to
converge asymptotically to local minima, and more importantly, at a better
rate as our algorithm possesses a better AMSE in comparison to the
aforementioned algorithms.

SO algorithms with non-convex objectives invariably suffer from the problem of
converging to stationary points that are not necessarily local minima. Ge et
al. (2015) suggests adding an additional noise term in the gradient estimate
when entering in the neighborhood of a stationary point that would ensure
escape from saddle points with high probability while working with an unbiased
SG estimate. In this work, however, we show that the inherent noise in the
biased stochastic gradient estimate helps one escape unstable equilibria such
as saddle points and local maxima without injecting additional noise. Another
approach tried out is to use a second-order method Zhu (2021), Balasubramanian
and Ghadimi (2018) to escape from unstable equilibria, but second-order
methods in general are more computationally expensive.

We now briefly review some prior work on non-asymptotic analyses in the
setting of a biased gradient oracle, which encapsulates the properties of a
gradient estimator based on SF. In Baes (2009), Devolder et al. (2014),
algorithms with a deterministic gradient oracle are presented. In Hu et al.
(2016), the authors analyze the rates achievable with inputs from a biased
stochastic gradient oracle, under a convexity assumption. In Bhavsar and
Prashanth (2022), the authors analyze SG algorithms for solving a non-convex
SO problem with inputs from a biased noisy gradient oracle. The rate that they
derive for a smooth objective matches the bound for our algorithm with a
balanced estimator. Further, unlike Bhavsar and Prashanth (2022), we derive an
improved non-asymptotic bound under a smoothness assumption for the noisy
observation $F$. This rate is $\mathcal{O}(1/\sqrt{N})$, where $N$ is the
number of iterations of our algorithm, and it matches the bound obtained for
GSF in Ghadimi and Lan (2013). In Balasubramanian and Ghadimi (2018), the
authors derive a non-asymptotic bound of $\mathcal{O}(1/\sqrt{N})$ for a
zeroth-order variant of the stochastic conditional gradient algorithm using
Gaussian perturbations. We derive a matching bound for our algorithm, which is
more efficient than the one in Balasubramanian and Ghadimi (2018) since their
algorithm requires solving an optimization problem in each iteration.

The rest of the paper is organized as follows: The framework for the
optimization problem and some preliminaries on GSF and TCSF are presented in
Section 2. Section 3 provides the main asymptotic and non-asymptotic
guarantees for TCSF while the convergence analysis is discussed in Section 4.
Section 5 presents simulation experiments that compare the performance of TCSF
with several algorithm.

##  2 SF gradient estimation with truncated Cauchy perturbations

We now present the idea behind a SF scheme for estimating the objective
function gradient proposed in Katkovnik and Kulchitsky (1972).
For a function $f:\mathbb{R}^{d}\to\mathbb{R}$, define the smoothed function
$g_{\delta}:\mathbb{R}^{d}\to\mathbb{R}$ as

$\displaystyle g_{\delta}(x)$
$\displaystyle=\mathbb{E}_{h_{\delta}(u)}[f(x+u)]=\mathbb{E}_{h_{\delta}(x-u)}[f(u)],\mbox{
}x\in\mathbb{R}^{d},$ (2)

where $h_{\delta}(x)$ is called the smoothing kernel or perturbation density
function. The parameter $\delta$ controls the degree of smoothness of
$g_{\delta}(x)$.
The SF scheme is useful especially if $f(x)$ is not well behaved, for
instance, if it has several stationary points in a narrow region. In such a
case, one may work with $g_{\delta}(x)$ as it would exhibit a smoother
behavior, and in general, it would be easier to compute it’s derivative as
opposed to $f(x)$. We have listed down the conditions from Rubinstein (1981)
which suggests that a density function needs to satisfy all of them to become
a smoothed function in the SF scheme.

#### Rubinstein’s Conditions for SF schemes

  1. (a)

$h_{\delta}:\mathbb{R}^{d}\to\mathbb{R}$ such that
$h_{\delta}(u)=\frac{1}{\delta^{d}}h(\frac{u}{\delta})$ is piecewise
differentiable with respect to $u$,

  2. (b)

$h_{\delta}(u)$ is a probability density function such that
$g_{\delta}(x)=\mathbb{E}_{h_{\delta}(u)}[f(x+u)]$ ,

  3. (c)

$\lim_{\delta\to 0}h_{\delta}(u)=\Delta(u)$, where $\Delta(u)$ is the Dirac-
Delta function,

  4. (d)

$\lim_{\delta\to 0}g_{\delta}(\cdot)=f(\cdot)$.

######  Definition 1.

[Truncated Cauchy Distribution] A r.v. is said to follow the truncated (to the
$\delta$-sphere) Cauchy distribution with mean vector zero and covariance
matrix $\Sigma=\delta^{2}\mathbb{I}_{d\times d}$ if its probability density
function has the following form (with $c_{1}$ being the normalization
constant):

$h_{\delta}(u)=\frac{\Gamma(\frac{d+1}{2})}{\pi^{\frac{d+1}{2}}c_{1}\delta^{d}}\frac{1}{(1+\frac{\lVert
u\rVert^{2}}{\delta^{2}})^{\frac{d+1}{2}}}\>\>\ \mbox{ for }\>\lVert
u\rVert^{2}\leq\delta^{2}.$ (3)

In the following remark we have shown how the truncated Cauchy distribution,
defined above, satisfies the aforementioned conditions.

######  Remark 1.

We now show that the truncated Cauchy distribution that we employ satisfies
the Rubinstein’s conditions (a)-(d) stated above. Hence, it is a valid
candidate to be used as a smoothing density functional in the SF algorithm.
Note that $h_{\delta}$ in our case is a truncated Cauchy distribution as
defined in 3 which is a probability density function. It is easy to see that
$h_{\delta}(u)=$ $\displaystyle\frac{1}{\delta^{d}}h(\frac{u}{\delta})$ where

$\displaystyle
h(\frac{u}{\delta})=\frac{\Gamma(\frac{d+1}{2})}{\pi^{\frac{d+1}{2}}c_{1}}\frac{1}{(1+\frac{\lVert
u\rVert^{2}}{\delta^{2}})^{\frac{d+1}{2}}}.$

Hence $h(u)$ denotes truncated Cauchy distribution with $\delta=1$. However
from 2 using this $h_{\delta}(u)$ one can write $g_{\delta}$ as expectation of
$f$ under $h_{\delta}(u)$ i.e.,
$g_{\delta}(x)=\mathbb{E}_{h_{\delta}(u)}[f(x+u)]$. We know that Dirac-Delta
function ($\Delta(u)$) is a measure whose value is $\infty$ at origin and $0$
otherwise with $\displaystyle\int_{-\infty}^{\infty}\Delta(u)du=1$. For
truncated Cauchy distribution, $\lim_{\delta\to 0}h_{\delta}(u)=\infty$ which
is a Delta function and $\lim_{\delta\to 0}g_{\delta}(\cdot)=f(\cdot)$ as
$\int_{\lVert u\rVert^{2}\leq 1}\Delta(u)du=1$. Thus, Rubinstein’s conditions
are satisfied in the case of the truncated Cauchy distribution.

A truncated Cauchy distribution has been considered in Chapter 6 of Bhatnagar
et al. (2013). However, the truncation there is not to the $\delta$-sphere,
and more importantly, unlike them, we derive asymptotic as well as non-
asymptotic rate results with the truncated Cauchy distribution specified
above.

One can intuitively interpret the effect of $\delta$ on smoothing as follows:
For smaller values of $\delta$, $g_{\delta}(x)$ is close to $f$. However, as
$\delta$ increases, $g_{\delta}(x)$ becomes smoother with fewer sharp
variations. As explained in Rubinstein (1981), the SF approach provides a
helpful way for approximating the gradient of any function $f$. In particular,
we have shown in the following proposition that for truncated Cauchy
smoothing, the derivative of $g_{\delta}(x)$ can be calculated by taking the
derivative of $h_{\delta}(x-u)$. This can be obtained via a simple application
of the Leibnitz rule and piece-wise differentiation property of
$h_{\delta}(u)$.

######  Proposition 1.

Let $g_{\delta}(x)$ be the smoothed function defined in 1 under truncated
Cauchy distribution. Then gradient of $g_{\delta}$ is

$\begin{split}\nabla
g_{\delta}(x)=\frac{1}{\delta}\mathbb{E}_{h(u)}\left[f(x+\delta
u)\frac{(d+1)u}{(1+\lVert u\rVert^{2})}\right].\\\ \end{split}$ (4)

###### Proof.

2 can be rewritten by considering $u=\delta v$ as follows:

$\begin{split}g_{\delta}(x)&=\frac{\Gamma(\frac{d+1}{2})}{\pi^{\frac{d+1}{2}}c_{1}}\int_{\left\lVert
v\right\rVert^{2}\leq 1}\frac{f(x+\delta v)}{(1+\left\lVert
v\right\rVert^{2})^{\frac{d+1}{2}}}dv.\\\ \end{split}$

Let’s rewrite the above equality by a change of variable $y=x+\delta v$ as
follows:

$\begin{split}g_{\delta}(x)&=\frac{\Gamma(\frac{d+1}{2})}{\delta^{d}\pi^{\frac{d+1}{2}}c_{1}}\int_{\left\lVert\frac{y-x}{\delta}\right\rVert^{2}\leq
1}\frac{f(y)}{(1+\left\lVert\frac{y-x}{\delta}\right\rVert^{2})^{\frac{d+1}{2}}}dy.\\\
\end{split}$

We apply the classic rule of differentiation over $x$ in the above equation,
and simplify as follows:

$\begin{split}\nabla g_{\delta}(x)&=\frac{1}{\delta}\int_{\lVert
u\rVert^{2}\leq 1}f(x+\delta u)\bigg{(}\frac{(d+1)u}{(1+\lVert
u\rVert^{2}}\bigg{)}h(u)du=\frac{1}{\delta}\mathbb{E}_{u}\bigg{[}\frac{f(x+\delta
u)(d+1)u}{(1+\lVert u\rVert^{2})}\bigg{]},\end{split}$

where $\mathbb{E}_{u}$ denotes expectation w.r.t. $h(u)$. ∎

In 4, $h(u)$ indicates that $h_{\delta}$ is computed at $\delta=1$. Notice
that the one-simulation smooth function under GSF scheme is specified in
Rubinstein (1981) by $\mathbb{E}_{h(u)}[f(x+\delta u)]$ and its derivative can
be approximated by ${\displaystyle\frac{1}{\delta}\mathbb{E}_{h(u)}[f(x+\delta
u)u]}$ where $u$ has the multivariate Gaussian distribution, a form also
studied in Bhatnagar and Borkar (2003) for the case when the objective
function is a long-run average cost. A finite-difference form of $\nabla
g_{\delta}(x)$ would be
${\displaystyle\frac{1}{\delta}\mathbb{E}_{h(u)}[(f(x+\delta u)-f(x))u]}$.

We now define the $\delta$-difference smoothed function
$f_{\delta}(x):\mathbb{R}^{d}\to\mathbb{R}$ for $f(x)$ (an approximation to
smoothed function) as below:

$\begin{split}f_{\delta}(x)&=\mathbb{E}_{h_{\delta}(u)}(f(x+u)-f(x))=\mathbb{E}_{h_{\delta}(x-u)}[f(u)-f(x-u)],\end{split}$
(5)

where $h_{\delta}(u)$ is smoothing kernel. The finite-difference gradient form
$\nabla g_{\delta}(x)$ for GSF , see Rubinstein (1981), is the same as $\nabla
f_{\delta}(x)$ due to the fact that $\mathbb{E}_{h(u)}(u)=0$ under Gaussian
smoothing kernel. As mentioned previously, we only have access to noisy
observations of the objective function. Thus, to solve the problem 1 using the
SF algorithm one can consider the gradient of $f_{\delta}(x)$ to be (with
$\xi$ having the same distribution as $\xi^{+}$) ${\displaystyle\nabla
f_{\delta}(x)=\frac{1}{\delta}\mathbb{E}_{u,\xi}[(F(x+\delta
u,\xi^{+})-F(x,\xi))u]}$.

It can similarly be shown as 4 that for $f\in\mathcal{C}_{L}^{1,1}$
(satisfying (A8)) and with the truncated Cauchy smooth kernel, the gradient of
$f_{\delta}$ can be expressed as follows:

$\begin{split}\nabla
f_{\delta}(x)=\frac{1}{\delta}\mathbb{E}_{h(u)}\left[(f(x+\delta
u)-f(x))\frac{(d+1)u}{(1+\lVert u\rVert^{2})}\right].\\\ \end{split}$ (6)

Algorithm 1 Truncated Cauchy Smoothed Functional (TCSF) Algorithm

Input: Initial point $x_{1}\in\mathbb{R}^{d}$, non-negative step-sizes
$\\{\gamma_{k}\\}_{k\geq 1}$, and smoothing parameter $\delta_{k}>0$.

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

Generate $u_{k}$ from 3, compute
$G(x_{k},\xi_{k}^{+},\xi_{k},u_{k},\delta_{k})$ using 8, and update
$x_{k+1}=x_{k}-\gamma_{k}G(x_{k},\xi_{k}^{+},\xi_{k},u_{k},\delta_{k}).$ (7)

end for

In this work, we propose the $\delta$-difference Cauchy smoothed functional
scheme. Now note that,
$\mathbb{E}_{h(u)}\left[\frac{u}{1+\lVert u\rVert^{2}}\right]\neq 0$ and so we
propose a two-simulation finite-difference gradient estimate instead of a one-
simulation estimate as it can be seen to have a lower bias.

For the case of noisy function observations, see (1), the SF gradient with
truncated Cauchy would simply be ${\displaystyle\nabla f_{\delta}(x)=}$
${\displaystyle\frac{1}{\delta}\mathbb{E}_{u,\xi^{+},\xi}\left[(F(x+\delta
u,\xi^{+})-F(x,\xi))\frac{(d+1)u}{(1+\lVert u\rVert^{2})}\right]}$ and a one-
sample gradient estimate would have the form ${\displaystyle\nabla
f(x)\approx}$ ${\displaystyle\bigg{(}\frac{F(x+\delta
u,\xi^{+})-F(x,\xi)}{\delta}\bigg{)}\frac{(d+1)u}{(1+\lVert u\rVert^{2})},}$
for $\delta>0$ small. When this estimate is used in a stochastic approximation
(SA) scheme, averaging would happen naturally. Thus, under Cauchy
perturbation, one can consider the estimate of $\nabla f(x_{k})$ for a given
parameter value $x_{k}$ (that in turn would be updated as per an SA scheme as
mentioned above) to be as follows:

$\displaystyle G(x_{k},\xi_{k}^{+},\xi_{k},u_{k},\delta_{k})$
$\displaystyle\stackrel{{\scriptstyle\triangle}}{{=}}{\displaystyle\bigg{(}\frac{F(x_{k}+\delta_{k}u_{k},\xi_{k}^{+})-F(x_{k},\xi_{k})}{\delta_{k}}\bigg{)}\frac{(d+1)u_{k}}{1+\lVert
u_{k}\rVert^{2}}}.$ (8)

In the above estimates, $\\{\xi_{k}\\}$ and $\\{\xi_{k}^{+}\\}$ constitute the
measurement noise r.vs in the two simulations and are assumed to be i.i.d.,
having a common distribution, and further independent of one another.
Algorithm 1 contains the details of the update procedure.

##### Motivation for truncation in Cauchy perturbations:

To get a bias bound for $G$ in (8), one needs the expectation of
$\frac{f(x+\delta_{k}u_{k})-f(x_{k})}{\delta_{k}}\frac{(d+1)u_{k}}{1+\lVert
u_{k}\rVert^{2}}+\frac{\eta_{k}^{+}-\eta_{k}}{\delta_{k}}\frac{(d+1)u_{k}}{1+\lVert
u_{k}\rVert^{2}}$, which can be obtained after a Taylor series expansion.
However, this expectation only exists if the mean and covariance matrix of the
perturbation distribution exist (see proof of Lemma 1 in Section 4.1). Thus we
incorporate truncation in the random perturbation.

##  3 Convergence Results

###  3.1 Asymptotic convergence

Let $\mathcal{F}_{k}=\sigma(x_{m},m\leq k;u_{m},\xi_{m}^{+},\xi_{m},m<k)$,
$k\geq 1$, denote a sequence of increasing sigma fields. Let
$\eta_{k}^{+}=F(x_{k}+\delta_{k}u_{k},\xi_{k}^{+})-f(x_{k}+\delta_{k}u_{k}),\;\;\eta_{k}=F(x_{k},\xi_{k})-f(x_{k})$,
respectively, where $\\{\delta_{k}\\}$ is a sequence of smoothing parameters
that diminishes to zero as $k\to\infty$. It is easy to see that
$\mathbb{E}[\eta_{k}^{+}-\eta_{k}|\mathcal{F}_{k}]=0$.

We now make the following assumptions as in Bhatnagar et al. (2013).

######  (A1).

The step-size $\gamma_{k}$ and the smoothing parameter $\delta_{k}$ are
positive for all $k$. Further, $\gamma_{k},\delta_{k}\rightarrow 0$ as
$k\to\infty$ and $\sum_{k}\gamma_{k}=\infty,\
\sum_{k}\left(\frac{\gamma_{k}}{\delta_{k}}\right)^{2}<\infty.$

######  (A2).

The function $f$ is three-times continuously differentiable with
$\rVert\nabla^{2}f(x)\lVert\leq B$ and $|\nabla^{3}_{i,j,k}f(x)|\leq B_{1}$
for $i,j,k=1,\ldots,N$ where $B,B_{1}\geq 0$

######  (A3).

There exist $\beta_{1},\beta_{2}>0$ such that $\forall k\geq 0$,
$\mathbb{E}|\eta_{k}|^{2}\leq\beta_{1},\mathbb{E}|\eta_{k}^{+}|^{2}\leq\beta_{1},\
\mathbb{E}|f(x_{k}\pm\delta_{k}u_{k})^{2}|\leq\beta_{2}$ and
$\mathbb{E}|f(x_{k})^{2}|\leq\beta_{2}$.

######  (A4).

$\sup_{k}\lVert x_{k}\rVert<\infty$ almost surely.

The above assumptions are commonly used in the convergence analysis of an SA
algorithm. In particular, (A1) is a standard SA requirement on the step-sizes.
(A2) ensures that the associated ODE is well-posed and helps in establishing
the asymptotic unbiasedness of the estimated gradient. The conditions in (A3)
ensure that the effect of noise can be ignored in asymptotic analysis of 7.
(A4) is a stability assumption that ensures convergence of 7 and is common to
the analysis of simultaneous perturbation based SG algorithm cf. Spall (1992),
Bhatnagar et al. (2013). If stability is hard to ensure, a common practice is
to project the iterate sequence onto a compact and convex set, and such a
scheme would fall under the realm of projected stochastic approximation.

The following lemma characterizes the relationship between the estimator $G$
and the true gradient of the objective $f$.

######  Lemma 1.

Under (A1)-(A4), we have almost surely

$\displaystyle\mathbb{E}[G(x_{k},\xi_{k}^{+},\xi_{k},u_{k},\delta_{k})|\mathcal{F}_{k}]=c_{2}\nabla
f(x_{k})+\delta_{k}w_{k},$ (9)

where $c_{2}=\mathbb{E}\left[\frac{(d+1)(u_{k}^{1})^{2}}{1+\lVert
u_{k}\rVert^{2}}\right]$, with $u_{k}^{1}$ denoting the first element of the
random vector $u_{k}$, and
$w_{k}=\mathbb{E}\left[\bigg{(}\frac{u_{k}^{T}\nabla^{2}f(\bar{x}^{+}_{k})u_{k}}{2}\bigg{)}\left.\frac{(d+1)u_{k}}{1+\lVert
u_{k}\rVert^{2}}\right|\mathcal{F}_{k}\right]$.

###### Proof.

See Section 4.1. ∎

Lemma 1 does not bound the bias in the gradient estimator directly. This is
because of the constant factor $c_{2}$ multiplying the gradient term on the
RHS of (9). However, the result in (9) is useful in establishing asymptotic
convergence of Algorithm 1 as it tracks the following ODE

$\dot{x(t)}=-c_{2}\nabla f(x(t)).$ (10)

In fact for $c_{2}>1$, we will obtain faster convergence, see 19 in Section
4.1 for a detailed argument.

######  Theorem 1.

(Strong Convergence): Assume (A1)-(A4) hold. The sequence $x_{k},k\geq 1$
governed by 7 satisfies

$x_{k}\to H\stackrel{{\scriptstyle\triangle}}{{=}}\\{x^{*}|\nabla
f(x^{*})=0\\}\textrm{ a.s. as }k\to\infty.$

###### Proof.

See Section 4.1.2. ∎

The following variant of (A3) is required to prove the asymptotic normality of
$x_{k}$.

######  (A5).

Assume (A3) holds. In addition, there exists $\sigma^{\prime}>0$ such that
$\mathbb{E}(\eta_{k}^{+}-\eta_{k})^{2}\to\sigma^{\prime 2}$ a.s. as
$k\to\infty$.

######  Theorem 2.

(Asymptotic Normality): Assume (A1)–(A5) hold. Also, assume that the set $H$
in Theorem 1 is a singleton $H=\\{x^{*}\\}$. Let
$\gamma_{k}=\gamma_{0}/k^{\alpha}$ and $\delta_{k}=\delta_{0}/k^{\phi}$, where
$\gamma_{0},\delta_{0}>0,\alpha\in(0,1]$ and $\phi\geq\alpha/6$. Furthermore,
let $\upsilon=\alpha-2\phi>0$ and $Q$ be the orthogonal matrix with
$Q\nabla^{2}f(x^{*})Q^{T}=\gamma_{0}^{-1}$diag$(\lambda_{1},\lambda_{2},\ldots,\lambda_{d})$,
with $\lambda_{1},\ldots,\lambda_{d}$ being the eigen-values of
$\nabla^{2}f(x^{*})$. Then,
$k^{\upsilon/2}(x_{k}-x^{*})\stackrel{{\scriptstyle
dist}}{{\to}}\mathcal{N}(\mu,QMQ^{T})$ as $k\to\infty$ a.s., where
$\mathcal{N}(\cdot,\cdot)$ denotes a multi-variate Gaussian distribution with
mean $\mu$ defined by

$\mu=\begin{cases}0&\text{if}\ 3\phi-\alpha/2>0,\\\
\bar{c}(\gamma_{0}\delta_{0}^{2}(\gamma_{0}\nabla^{2}f(x^{*})-\frac{1}{2}\upsilon^{+}\mathbb{I})^{-1}T)&\text{if}\
3\phi-\alpha/2=0.\end{cases}$

In the above, $\mathbb{I}$ is identity matrix,$\bar{c}=\mathbb{E}[\lVert
u\rVert^{4}]$ , $\upsilon^{+}=\upsilon$ if $\alpha=1$ and $0$ if $\alpha<1,\
T$ is a d-dimensional vector whose $i$-th element given by
$-\frac{1}{6}[\nabla^{3}_{iii}f(x^{*})+3\sum_{j=1,j\neq
i}\nabla^{3}_{jji}f(x^{*})]$ and the matrix is defined as
$M=\frac{\gamma^{2}_{0}\sigma^{\prime
2}}{4\delta^{2}_{0}}diag((2\lambda_{1}-\upsilon^{+})^{-1},\ldots,(2\lambda_{d}-\upsilon^{+})^{-1}).$

###### Proof.

See Section 4.1.3. ∎

######  Remark 2.

For the case of general $H$ (not necessarily singleton as in the statement of
Theorem 2, the result will continue to hold for the particular $x^{*}$ in
whose neighborhood the parameter lies after a sufficiently large number of
iterations.

From Theorem 2 we can say that $k^{\upsilon/2}(x_{k}-x^{*})$ is asymptotically
Gaussian for TCSF algorithm. In addition, the maximum possible value of
$\upsilon=2/3$ can be obtained by fixing $\alpha=1$ and $\phi=1/6$ in
$\upsilon=\alpha-2\phi$. Fixing $\upsilon=2/3$, we obtain the best possible
asymptotic convergence rate of $k^{-1/3}$. We define the asymptotic mean
square error (AMSE), cf. Spall (1992), by
$AMSE(\gamma_{0},\delta_{0})=\mu^{T}\mu+trace(QMQ^{T})$, where
$\gamma_{0},\delta_{0}$ are constant step-size and smoothing constant
respectively. It can be shown under the condition given in Gerencser (1999)
that $AMSE(\gamma_{0},\delta_{0})$ coincides with
$k^{\upsilon}\mathbb{E}\lVert x_{k}-x^{*}\rVert^{2}$.

We now compare the AMSE of our algorithm with two other well-known random
perturbation gradient estimation schemes, namely GSF and SPSA

To make the comparison fair, we follow the approach of Chin (1997) and set
$\gamma_{k}$ and $\delta_{k}$ uniformly. In particular, we use step-size
$\gamma_{k}=\gamma_{0}/k$, where $\gamma_{0}\geq\upsilon/2\lambda_{0}$, with
$\lambda_{0}$ denoting the minimum eigenvalue of $\nabla^{2}f(x^{*})$.
Further, we set $\delta_{k}=\frac{\delta_{0}}{k^{1/6}}$. It can then be shown
that

$AMSE(\gamma_{0},\delta_{0})=\Big{(}\bar{c}\delta^{2}_{0}\gamma_{0}\lVert\Phi
T\rVert\Big{)}^{2}+\frac{1}{\delta^{2}_{0}}trace(\Phi P),$ (11)

where $T$ is as defined in Theorem 2,
$\Phi=(\gamma_{0}\nabla^{2}f(x^{*})-\frac{1}{2}\upsilon^{+}\mathbb{I})^{-1}$
and $P=\frac{\sigma^{\prime 2}}{4}\mathbb{I}$.

######  Remark 3 (Comparing with GSF).

For $u_{k}$ following $\mathcal{N}(0,\mathbb{I})$, we have
$\mathbb{E}(u_{k}^{i})^{4}=3$. Hence, the ratio of AMSE of GSF with that of
TCSF algorithm is given by

$\displaystyle\frac{AMSE_{GSF}(\gamma_{0},\delta_{0})}{AMSE_{TCSF}(\gamma_{0},\delta_{0})}$
$\displaystyle=\frac{\Big{(}3\delta^{2}_{0}\gamma_{0}\lVert\Phi
T\rVert\Big{)}^{2}+\frac{1}{\delta^{2}_{0}}trace(\Phi
P)}{\Big{(}\bar{c}\delta^{2}_{0}\gamma_{0}\lVert\Phi
T\rVert\Big{)}^{2}+\frac{1}{\delta^{2}_{0}}trace(\Phi P)}.$

where $\Phi,\ P$ $T$ are specified in (11). Recall that
$\bar{c}=\mathbb{E}[\lVert u\rVert^{4}]$, and $u$ is restricted to the unit
sphere, implying $\bar{c}\leq 1$. Thus, the AMSE of TCSF is clearly better
than that of GSF.

######  Remark 4 (Comparing with SPSA).

Comparing the AMSE of TCSF with SPSA, we obtain

$\displaystyle\frac{AMSE_{SPSA}(\gamma_{0},\delta_{0})}{AMSE_{TCSF}(\gamma_{0},\delta_{0})}$
$\displaystyle=\frac{\Big{(}\delta^{2}_{0}\gamma_{0}\lVert\Phi
T\rVert\Big{)}^{2}+\frac{1}{\delta^{2}_{0}}trace(\Phi
P)}{\Big{(}\bar{c}\delta^{2}_{0}\gamma_{0}\lVert\Phi
T\rVert\Big{)}^{2}+\frac{1}{\delta^{2}_{0}}trace(\Phi P)}\geq 1.$

The last inequality follows from the fact that $\bar{c}\leq 1$. Thus, the AMSE
of TCSF is at least as good as that of SPSA, and would be better if
$\bar{c}<1$. For the case of $d=1$, $\bar{c}<1$, as shown in Staneski (1990).
On the other hand, it is difficult to obtain a closed-form expression for
$\bar{c}$ when $d>1$.

From the analysis above, it is apparent that, from an asymptotic convergence
rate viewpoint, our algorithm outperforms GSF and SPSA, which are two popular
gradient estimation schemes.
So far, we have provided theoretical guarantees that establish convergence to
a stationary point of the objective function $f$. However, this result is not
sufficient in a non-convex optimization setting since local maxima and saddle
points are also stationary points in addition to local minima. The aim here is
to converge to a local minimum, or avoid traps such as saddle points/local
maxima. A usual trick to achieve this objective is to add extraneous noise, so
that the algorithm does not converge to an unstable equilibrium. We now
establish that our algorithm naturally avoids traps owing to the noise in the
gradient estimator. For this result, we require the following additional
assumption.

######  (A6).

Assume the condition in (A5) holds and $c_{9}$ is such that
$\mathbb{E}|\eta_{k}^{+}-\eta_{k}|\geq c_{9}$.

The assumption above ensures that the noise is rich in all directions. Note
that, one can rewrite the update rule 7 as

$x_{k+1}=x_{k}-\gamma_{k}G(x_{k},u_{k},\delta_{k})-\epsilon_{k},$ (12)

with ${\displaystyle
G(x_{k},u_{k},\delta_{k})=\frac{f(x_{k}+\delta_{k}u_{k})-f(x_{k})}{\delta_{k}}\frac{(d+1)u_{k}}{1+\lVert
u_{k}\rVert^{2}}}$ and
${\displaystyle\epsilon_{k}=\gamma_{k}\frac{\eta_{k}^{+}-\eta_{k}}{\delta_{k}}\frac{(d+1)u_{k}}{1+\lVert
u_{k}\rVert^{2}}}$.
Here $\mathbb{E}[\epsilon_{k}|\mathcal{F}_{k}]=0$. Note that 12 is equivalent
to 23 in Section 4.1 by considering $-\epsilon_{k}=\mu_{k}$ and $-G=Y$. One
can ensure avoidance of traps if the increment of $\epsilon_{k}$ in any
direction is of order $1/n^{\gamma}$, i.e.,
$\displaystyle\mathbb{E}\left[(\epsilon_{k}\cdot\theta)^{+}|\mathcal{F}_{k}\right]\geq
c_{7}/k^{\gamma}$ (Theorem 7) for every unit vector $\theta$. We establish
that our algorithm satisfies
$\displaystyle\mathbb{E}\left[(\epsilon_{k}\cdot\theta)^{+}|\mathcal{F}_{k}\right]\geq
c_{9}c_{10}/2$ with $c_{9}=\mathcal{O}(1/k^{\gamma})$ (see (Section 4.1.4) for
details).

######  Proposition 2.

Under (A1)-(A6), $x_{k}$ generated by Algorithm 1, converges to a local
minimum a.s.

###### Proof.

See Section 4.1.4. ∎

######  Remark 5.

From Proposition 2 we can justify that 7 avoids saddle points and local
maxima. To the best of our knowledge, a similar result is not available for
the GSF algorithm. The latter algorithm has been shown to converge to a
stationary point in Bhatnagar et al. (2013), and a non-asymptotic convergence
rate for the same is available in Ghadimi and Lan (2013). In Ge et al. (2015),
Jin et al. (2017), the authors suggest adding extraneous noise to avoid traps
for a SG algorithm. In contrast, we show that a noisy gradient estimation
scheme would naturally avoid traps, obviating the need for extraneous noise
addition.

###  3.2 Non-asymptotic convergence

The non-asymptotic analysis below establishes convergence to an approximate
stationary point as with Ghadimi and Lan (2013), Bhavsar and Prashanth (2022).

######  Definition 2.

For a non-convex function $f$, $\bar{x}$ is said to be an
$\epsilon$-stationary point to the problem 1 if it satisfies
$\mathbb{E}\left[\left\lVert\nabla
f(\bar{x})\right\rVert^{2}\right]\leq\epsilon$.

For non-asymptotic analysis we make the following assumption.

######  (A7).

Let $\sigma^{2}>0$ such that $\mathbb{E}_{\xi}[\lVert\nabla F(x,\xi)-\nabla
f(x)\rVert^{2}]\leq\sigma^{2}$, for all $x\in\mathbb{R}^{d}$.

######  (A8).

The function $f:\mathbb{R}^{d}\to\mathbb{R}$ is differentiable and the
gradient of $f$ is Lipschitz continuous, i.e., $\left\lVert\nabla f(x)-\nabla
f(y)\right\rVert\leq L\left\lVert x-y\right\rVert$, $\forall
x,y\in\mathbb{R}^{d}$, where $L>0$ is the Lipschitz constant.

We now present a non-asymptotic bound for a randomized variant of TCSF
algorithm in the spirit of Ghadimi and Lan (2013).

######  Theorem 3.

Assume (A3), (A7) and (A8). Suppose Algorithm 1 has the following as the step-
sizes and smoothing parameters:

$\gamma_{k}\stackrel{{\scriptstyle\triangle}}{{=}}\min\left\\{\dfrac{c_{2}}{L},\dfrac{1}{N^{2/3}}\right\\}\
,\
\delta_{k}\stackrel{{\scriptstyle\triangle}}{{=}}\delta=\dfrac{1}{N^{1/6}},k=1,\ldots,N,$

where $c_{2}$ is as specified in Lemma 1. Let $x_{R}$ be picked uniformly at
random from $\\{x_{1},\ldots,x_{N}\\}$. Then

$\begin{split}\mathbb{E}\left[\left\lVert\nabla
f(x_{R})\right\rVert^{2}\right]&\leq
2D\Big{(}\frac{L}{c_{2}N}+\frac{1}{N^{1/3}}\Big{)}+\frac{2Bc_{2}^{\prime}}{c_{2}N^{1/6}}+\frac{L}{2c_{2}}\left(\frac{2Bc_{2}^{\prime}}{N^{5/6}}+\frac{dc_{2}^{\prime
2}}{N}+\frac{C^{\prime\prime}}{N^{1/3}}\right),\end{split}$

where $C^{\prime\prime}=\mathbb{E}\left[\frac{\lVert
u_{k}\rVert^{2}(d+1)^{2}}{(1+\lVert
u_{k}\rVert^{2})^{2}}\right]\big{[}2(\beta_{1}+\beta_{2})\big{]}$ ,
$c_{2}^{\prime}=(d+1)B$ and $B$ is as specified in (A2).

###### Proof.

See Section 4.2.1. ∎

From Theorem 3, we can conclude that the non-asymptotic bound is
$\mathcal{O}(N^{-1/6})$, that is weaker than the asymptotic convergence rate
of $\mathcal{O}(N^{-1/3})$ obtained in the previous section. This gap can
however be fixed by using a balanced estimator, i.e.

$\begin{split}\tilde{G}(x_{k},\xi_{k}^{+},\xi_{k},u_{k},\delta_{k})&\stackrel{{\scriptstyle\triangle}}{{=}}{\displaystyle\bigg{(}\frac{F(x_{k}+\delta_{k}u_{k},\xi_{k}^{+})-F(x_{k}-\delta_{k}u_{k},\xi_{k}^{-})}{2\delta_{k}}\bigg{)}\frac{(d+1)u_{k}}{1+\lVert
u_{k}\rVert^{2}}}.\end{split}$ (13)

In this case, we obtain the following bound (see Lemma 7 in Section 4.3 for
details):

$\mathbb{E}[\tilde{G}(x_{k},\xi_{k}^{+},\xi_{k},u_{k},\delta_{k})|\mathcal{F}_{k}]=c_{2}\nabla
f(x_{k})+\mathcal{O}(\delta^{2}_{k})\mathbf{1}_{d}.$

Notice that the bound is $\mathcal{O}(\delta_{k}^{2})$ instead of
$\mathcal{O}(\delta_{k})$ in the one-sided estimator (8), since the second-
order terms cancel in the balanced estimator. With a
$\mathcal{O}(\delta_{k}^{2})$ bias bound, we obtain in the following theorem a
non-asymptotic bound of $\mathcal{O}(N^{-1/3})$ using the balanced estimator.
From this bound it is observed that $\mathcal{O}(1/\epsilon^{3})$ number of
iterations of the TCSF algorithm guarantee convergence to an
$\epsilon$-stationary point of 1.

######  Theorem 4.

Assume (A3), (A7) and (A8). Suppose Algorithm 1 is running with step-sizes
$\gamma_{k}$ and $\tilde{G}(x_{k},\xi_{k}^{+},\xi_{k},u_{k},\delta_{k})$
instead of $G(x_{k},\xi_{k}^{+},\xi_{k},u_{k},\delta_{k})$. The smoothing
parameter $\delta_{k}=\delta,\forall k$, and step-sizes $\gamma_{k}$ is chosen
as defined in Theorem 3. Let $x_{R}$ denote a point picked uniformly at random
from $\\{x_{1},\ldots,x_{N}\\}$. Then under probability distribution 26

$\begin{split}\mathbb{E}\left[\left\lVert\nabla
f(x_{R})\right\rVert^{2}\right]&\leq\Big{(}\frac{2DL}{c_{2}N}+\frac{2D}{N^{1/3}}\Big{)}+\frac{2Bc_{2}^{\prime\prime}}{c_{2}N^{1/3}}+\frac{L}{2c_{2}}\left(\frac{2Bc_{2}^{\prime\prime}}{N}+\frac{dc_{2}^{\prime\prime
2}}{N^{4/3}}+\frac{C^{\prime\prime}}{N^{1/3}}\right),\end{split}$

where $C^{\prime\prime}$ is as defined in Section 4.1, $c_{2}^{\prime\prime}$
is defined Lemma 7 and $B$ is the same as in (A2).

###### Proof.

See Section 4.3.1 ∎

For the non-asymptotic bounds presented above, we assumed that the objective
function is smooth. Instead, if we make the stronger assumption that the noisy
observation $F$ is smooth, then we can obtain a better non-asymptotic bound.
We make this claim precise in the following.

######  (A9).

The function $F$ is Lipschitz continuous in the first argument, uniformly
w.r.t the second, i.e., for any given $\xi$, $\left\lVert\nabla
F(x,\xi)-\nabla F(y,\xi)\right\rVert\leq L\left\lVert x-y\right\rVert$,
$\forall x,y\in\mathbb{R}^{d}$ almost surely.

Assumptions (A7) and (A9) imply (A8), since

$\left\lVert\nabla f(x)-\nabla
f(y)\right\rVert\leq\mathbb{E}_{\xi}\left\lVert\nabla F(x,\xi)-\nabla
F(y,\xi)\right\rVert\leq L\left\lVert x-y\right\rVert,\;\;\forall
x,y\in\mathbb{R}^{d}.$

We consider now the following variant of the gradient estimator in (8):

$\displaystyle G(x_{k},\xi_{k},u_{k},\delta_{k})$
$\displaystyle\stackrel{{\scriptstyle\triangle}}{{=}}{\displaystyle\bigg{(}\frac{F(x_{k}+\delta_{k}u_{k},\xi_{k})-F(x_{k},\xi_{k})}{\delta_{k}}\bigg{)}\frac{(d+1)u_{k}}{1+\lVert
u_{k}\rVert^{2}}}.$ (14)

Notice that both function observations in this case use the same noise factor
$\xi_{k}$. Such a setting is possible when noise is added using common random
numbers, for instance, in computer simulations. In this setting,
$G(x_{k},\xi_{k},u_{k},\delta_{k})$ defined in (14) satisfies

$\mathbb{E}_{\xi,h(u)}\left[G(x_{k},\xi_{k},u_{k},\delta)\right]=\mathbb{E}_{h(u)}\left[\mathbb{E}_{\xi}\left[G(x_{k},\xi_{k},u_{k},\delta)\right]\right]=\nabla
f_{\delta}(x_{k}).$

We now provide a non-asymptotic bound of the order $\mathcal{O}(N^{-1/2})$ for
Algorithm 1 under the additional assumptions listed above.

######  Theorem 5.

Assume (A7), (A9). Suppose Algorithm 1 runs with step-sizes $\gamma_{k}$ and
the smoothing parameter $\delta_{k}=\delta,\forall k$, chosen as follows:

$\gamma_{k}=\min\left\\{\frac{1}{2Lc_{13}},\frac{1}{c_{13}\sigma\sqrt{N}}\right\\},\
\delta=\frac{1}{L\sqrt{dNc_{13}}},$

where $c_{13}=\frac{4c_{11}c_{12}}{d+1}$, with $c_{11}$ defined in Theorem 2
and $c_{12}$ is the Frobenius norm of the generalized inverse of the matrix
$\mathbb{E}_{u}\left(\frac{(d+1)uu^{T}}{1+\lVert u\rVert^{2}}\right)$. Let
$x_{R}$ denote a point picked uniformly at random from
$\\{x_{1},\ldots,x_{N}\\}$. Then

$\mathbb{E}\left[\|\nabla
f(x_{R})\|^{2}\right]\leq\frac{c_{14}}{N}+\frac{2\sigma L}{\sqrt{N}},\mbox{
where }c_{14}=2L+c_{13}\sigma+8LBc_{13}^{3/2}+\frac{\sigma
Bc_{13}^{3/2}B}{L}+\frac{2c_{13}d+1}{L}.$

###### Proof.

See Section 4.2. ∎

######  Remark 6.

In Theorem 5, we notice that the convergence rate is $\mathcal{O}(1/\sqrt{N})$
or equivalently $\mathcal{O}(1/\epsilon^{2})$ number of iterations are needed
to find an $\epsilon$-stationary point of 1. This rate is better than the one
obtained in Theorem 3 and this is a consequence of (A9) which ensures $F$ is
smooth. In Ghadimi and Lan (2013), for GSF, the number of iterations to find
an $\epsilon$-stationary point is bounded by $\mathcal{O}(1/\epsilon^{2})$,
and our bound matches their result. The advantage with our algorithm is that
it outperforms the GSF algorithm empirically. In the next section, we provide
some examples to validate this claim.

We will provide in the following theorem a non-asymptotic bound of order
$\mathcal{O}(N^{-}{1/2})$ for Algorithm 1 with balanced estimator by assuming
sample performance is smooth.

######  Theorem 6.

Assume (A2), (A7) and (A9). Suppose Algorithm 1 is running with
$\tilde{G}(x_{k},\xi_{k}^{+},\xi_{k},u_{k},\delta_{k})$ instead of
$G(x_{k},\xi_{k}^{+},\xi_{k},u_{k},\delta_{k})$. Let the smoothing parameter
$\delta_{k}=\delta,\forall k$ and step-sizes $\gamma_{k}$ is chosen as
follows:

$\begin{split}\gamma_{k}=\min\left\\{\frac{c_{2}}{2c^{2}_{11}L},\frac{1}{N^{1/2}}\right\\},\
\delta=\frac{1}{N^{1/2}}.\end{split}$ (15)

Here $c_{2}$ is same as in Lemma 1. Let $x_{R}$ denote a point picked
uniformly at random from $\\{x_{1},\ldots,x_{N}\\}$. Then under probability
distribution $\mathbb{P}_{R}(k)=\frac{1}{N}$

$\begin{split}\mathbb{E}\left[\left\lVert\nabla
f(x_{R})\right\rVert^{2}\right]&\leq\Big{(}\frac{2DL}{Nc_{2}^{2}}+\frac{2D}{N^{1/2}}\Big{)}+\frac{2Bc_{2}^{\prime\prime}}{c_{2}N}+\frac{L}{2c_{2}N^{1/2}}\left(4c_{11}\sigma^{2}+\frac{2d^{2}c_{2}^{\prime\prime
2}}{N^{2}}\right),\end{split}$

where $c_{2}^{\prime\prime}$ is defined Lemma 7 and $B$ is the same as in
(A2).

###### Proof.

See Section 4.3.2 ∎

##  4 Convergence proofs

###  4.1 Asymptotic convergence proofs

####  4.1.1 Proof of Lemma 1

In the lemma below, we state and prove a bound on the moments of a truncated
Cauchy r.v.

######  Lemma 2.

Let u be a truncated Cauchy r.v., then for any positive integer r, we have

$\mathbb{E}_{u}[\lVert u\rVert^{2r}]\leq\frac{c_{11}}{(r+d)},$ (16)

where $c_{11}=\frac{2\Gamma(\frac{d+1}{2})}{\sqrt{\pi}\Gamma(d/2)c_{1}}.$

###### Proof.

$\begin{split}\mathbb{E}_{u}[\lVert
u\rVert^{2r}]&=\frac{\Gamma(\frac{d+1}{2})}{\pi^{\frac{d+1}{2}}c_{1}}\int_{\lVert
u\rVert^{2}\leq 1}\lVert u\rVert^{2r}\frac{1}{(1+\lVert
u\rVert^{2})^{\frac{d+1}{2}}}du\\\
&=\frac{\Gamma(\frac{d+1}{2})}{\pi^{\frac{d+1}{2}}c_{1}}\int_{0}^{1}\int_{\lVert
u\rVert^{2}=R}R^{r}\frac{1}{(1+R)^{\frac{d+1}{2}}}dudR\\\
&=\frac{\Gamma(\frac{d+1}{2})}{\pi^{\frac{d+1}{2}}c_{1}}\int_{0}^{1}R^{r}\frac{1}{(1+R)^{\frac{d+1}{2}}}\frac{2\pi^{d/2}}{\Gamma(d/2)}R^{d-1}dR\\\
&=c_{11}\int_{0}^{1}\frac{R^{r+d-1}}{(1+R)^{\frac{d+1}{2}}}dR\\\ &\leq
c_{11}\int_{0}^{1}R^{r+d-1}dR\\\ &=\frac{c_{11}}{(r+d)}.\end{split}$

∎

#### Proof of Lemma 1

###### Proof.

Using Taylor series expansions, we obtain

$\begin{split}f(x_{k}+\delta_{k}u_{k})&=f(x_{k})+\delta_{k}u_{k}^{T}\nabla
f(x_{k})+\frac{\delta_{k}^{2}}{2}u_{k}^{T}\nabla^{2}f(\bar{x}^{+}_{k})u_{k},\end{split}$
(17)

where $\bar{x}^{+}_{k}$ is on the line segment between $x_{k}$ and
$x_{k}+\delta_{k}u_{k}$. Using (17), we have

$\displaystyle G(x_{k},\xi_{k}^{+},\xi_{k},u_{k},\delta_{k})$
$\displaystyle=\frac{f(x+\delta_{k}u_{k})-f(x_{k})}{\delta_{k}}\frac{(d+1)u_{k}}{1+\lVert
u_{k}\rVert^{2}}+\frac{(\eta_{k}^{+}-\eta_{k})}{\delta_{k}}\frac{(d+1)u_{k}}{1+\lVert
u_{k}\rVert^{2}},$ $\displaystyle=\frac{(d+1)u_{k}u_{k}^{T}\nabla
f(x_{k})}{1+\lVert
u_{k}\rVert^{2}}+\bigg{(}\frac{u_{k}^{T}\nabla^{2}f(\bar{x}^{+}_{k})u_{k}}{2}\bigg{)}\frac{(d+1)u_{k}}{1+\lVert
u_{k}\rVert^{2}}\delta_{k}+\frac{(\eta_{k}^{+}-\eta_{k})}{\delta_{k}}\frac{(d+1)u_{k}}{1+\lVert
u_{k}\rVert^{2}},$

Let $u_{k}^{i}$ denote the $i$-th element of the random vector $u_{k}$. Then,
the off-diagonal elements of
$\displaystyle\mathbb{E}_{u}\bigg{[}\frac{(d+1)u_{k}u_{k}^{T}}{1+\lVert
u_{k}\rVert^{2}}\bigg{]}$ satisfy
$\displaystyle\mathbb{E}_{u}\bigg{[}\frac{(d+1)u_{k}^{i}u_{k}^{j}}{1+\lVert
u_{k}\rVert^{2}}\bigg{]}=0,$ since $\displaystyle\frac{1}{1+\lVert
u_{k}\rVert^{2}}\in(0,1]$ is upper bounded by a fixed quantity, and for $i\neq
j$, $\mathbb{E}_{u}[u_{k}^{i}u_{k}^{j}]=0$. Hence

$\mathbb{E}_{u}\left[\frac{(d+1)u_{k}u_{k}^{T}}{1+\lVert
u_{k}\rVert^{2}}\right]=c_{2}\mathbb{I},$ (18)

where $c_{2}$ is as defined in the lemma statement.

Notice that $\mathbb{E}_{u}\left[(\eta_{k}^{+}-\eta_{k})\frac{u_{k}}{1+\lVert
u_{k}\rVert^{2}}\middle|\ \mathcal{F}_{k}\right]=0$, since
$E[\eta_{k}^{+}-\eta_{k}\mid\mathcal{F}_{k}]=0$, and $u_{k}$ is independent of
$\mathcal{F}_{k}$. Thus we obtain

$\displaystyle\mathbb{E}[G(x_{k},\xi_{k}^{+},\xi_{k},u_{k},\delta_{k})|\mathcal{F}_{k}]=c_{2}\nabla
f(x_{k})+\delta_{k}w_{k}.$

∎

####  4.1.2 Proof of Theorem 1

For the proof of Theorem 1, we require the notion of Lyapunov stability, which
we define next.

######  Definition 3.

A continuously differentiable function $V\colon\mathbb{R}^{d}\to[0,\infty)$ is
said to be a Lyapunov function for an ODE $\dot{y}=f(y)$ with set of
equilibrium points $\mathit{H}$ if it satisfies the properties below.

  1. 1.

$\lim_{\lVert x\rVert\to+\infty}V(x)=\infty$.

  2. 2.

The inner product of $f(y)$ with $V(y)$ can take the following values:

$\langle f(y),\nabla V(y)\rangle=\begin{cases}0&\text{if}\ y\in\mathit{H},\\\
<0&\text{otherwise}.\\\ \end{cases}$

######  Lemma 3.

Consider the function in 1 such that $f(\cdot)\geq c^{\prime}$ where
$c^{\prime}$ is a negative real number and $H$ be the set of equilibrium
points of the ODE $\dot{x}(t)=-\nabla f(x(t))$, i.e., $H=\\{x(t):\nabla
f(x(t))=0\\}$ where $f(x)$ is defined as in 1. Then $x(t)\to H$ as
$t\rightarrow\infty$.

###### Proof.

Let $g(\cdot)=f(\cdot)-c^{\prime}\geq 0$. Further

$\begin{split}\frac{dg(x(t))}{dt}&=\langle\nabla f(x(t)),\dot{x}(t)\rangle\\\
&=\langle\nabla f(x(t)),-\nabla f(x(t))\rangle\\\ &=-\lVert\nabla
f(x(t))\rVert^{2}.\end{split}$

Thus, ${\displaystyle\frac{dg(x(t))}{dt}<0}$ for $x(t)\not\in H$ and is 0
otherwise. Thus, $g$ serves as a Lyapunov function for the above (gradient)
ODE. The claim follows. ∎

#### Proof of Theorem 1

###### Proof.

From 7 and Lemma 1, we have

$\begin{split}x_{k+1}&=x_{k}-\gamma_{k}(\mathbb{E}_{u}[G(x_{k},\xi_{k}^{+},\xi_{k},u_{k},\delta_{k})|\mathcal{F}_{k}]+M_{k}),\\\
&=x_{k}-\gamma_{k}(c_{2}\nabla
f(x_{k})+c_{2}^{\prime}\delta_{k}\mathbf{1}_{d}+M_{k}).\\\ \end{split}$ (19)

where
$M_{k}=G(x_{k},\xi_{k}^{+},\xi_{k},u_{k},\delta_{k})-\mathbb{E}_{u}[G(x_{k},\xi_{k},\xi_{k}^{+},u_{k},\delta)|\mathcal{F}_{k}]$,
$k\geq 0$, is a martingale difference sequence. Further, $\gamma_{k},k\geq 1$
satisfies (A1). The update rule 19 thus tracks the ODE 10. However the ODE 10
has the same equilibria as the ODE

$\dot{x}=-\nabla f(x).$

In fact, if the constant $c_{2}>1$, the ODE will have a faster speed of
convergence to $x^{*}$. Note that $c_{2}^{\prime}\delta_{k}\mathbf{1}_{d}$ is
equivalent to the bias vector. Here each term of
$\mathcal{O}(\delta_{k})\mathbf{1}_{d}$ goes to zero as $k\to\infty$. From
Kushner and Clark (1978a) we can directly conclude the convergence of the
above algorithm using the assumptions (a)-(d) below that are taken from
Kushner and Clark (1978a).

  1. (a)

$\nabla f(x)$ is a Lipschitz continuous function.

  2. (b)

The bias sequence $\delta_{k},k\geq 1$ is almost surely bounded and with
$\delta_{k}\to 0$ almost surely as $k\to\infty$.

  3. (c)

The step-sizes $\gamma_{k},k\geq 1$, satisfy $\gamma_{k}\to 0$ as $k\to\infty$
and $\sum_{k}\gamma_{k}=\infty$.

  4. (d)

The Martingale sequence $M_{k}$ satisfies the following: $\forall\zeta>0$,

$\lim_{k\to\infty}\mathbb{P}\left(\sup_{m\geq
k}\lVert\sum_{i=k}^{m}\gamma_{i}M_{i}\rVert\geq\zeta\right)=0.$

Assumptions $(a),(b)$ and $(c)$ above directly follow from (A1) and (A2). We
now verify Assumption $(d)$ above. Recall the Doob’s martingale inequality,
i.e.

$\mathbb{P}\left(\sup_{m\geq 0}\lVert
Z_{m}\rVert\geq\zeta\right)\leq\frac{1}{\zeta^{2}}\lim_{m\to\infty}\mathbb{E}\lVert
Z_{m}\rVert^{2}.$

By considering $Z_{m}=\sum_{i=0}^{m-1}\gamma_{i}M_{i}$, we obtain

$\displaystyle\mathbb{P}\left(\sup_{m\geq
k}\lVert\sum_{i=k}^{m}\gamma_{i}M_{i}\rVert\geq\zeta\right)$
$\displaystyle\leq\frac{1}{\zeta^{2}}\lim_{m\to\infty}\mathbb{E}\lVert\sum_{i=k}^{m}\gamma_{i}M_{i}\rVert^{2}$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}\frac{1}{\zeta^{2}}\sum_{i=k}^{\infty}\gamma_{k}^{2}\mathbb{E}\lVert
M_{i}\rVert^{2},$ (20)

where the first inequality in $(\ref{def:eq6})$ follows from the fact that
$\mathbb{E}[M_{i}M_{j}]=\mathbb{E}[M_{i}\mathbb{E}[M_{j}]|\mathcal{F}_{j}]=0$.
Now using the identity $\mathbb{E}\lVert
X-\mathbb{E}[X|\mathcal{F}]\rVert^{2}\leq\mathbb{E}[\lVert X\rVert^{2}]$ for a
r.v. $X$ one can rewrite $\mathbb{E}\lVert M_{k}\rVert^{2}\leq\mathbb{E}\lVert
G(x_{k},\xi_{k}^{+},\xi_{k},u_{k},\delta_{k})\rVert^{2}$. Hence

$\displaystyle\mathbb{E}\lVert M_{k}\rVert^{2}$
$\displaystyle\leq\mathbb{E}\lVert
G(x_{k},\xi_{k}^{+},\xi_{k},u_{k},\delta_{k})\rVert^{2}$
$\displaystyle=\mathbb{E}\left[\frac{\lVert
u_{k}\rVert^{2}(d+1)^{2}}{(1+\lVert
u_{k}\rVert^{2})^{2}}\left(\frac{\eta_{k}^{+}-\eta_{k}}{\delta_{k}}\right)^{2}\right]+\mathbb{E}\left[\frac{\lVert
u_{k}\rVert^{2}(d+1)^{2}}{(1+\lVert
u_{k}\rVert^{2})^{2}}\left(\frac{f(x_{k}+\delta_{k}u_{k})-f(x_{k})}{\delta_{k}}\right)^{2}\right]$
$\displaystyle\leq\frac{C^{\prime\prime}}{\delta_{k}^{2}}.$ (21)

Here $C^{\prime\prime}=\mathbb{E}\left[\frac{\lVert
u_{k}\rVert^{2}(d+1)^{2}}{(1+\lVert
u_{k}\rVert^{2})^{2}}\right]\big{[}2(\beta_{1}+\beta_{2})\big{]}<\infty$. The
last inequality follows from (A3). Moreover, from construction, the truncated
Cauchy distribution has finite moments. Plugging the above inequality in
Section 4.1,we obtain

$\begin{split}\lim_{k\to\infty}\mathbb{P}\left(\sup_{m\geq
k}\lVert\sum_{i=k}^{m}\gamma_{i}M_{i}\rVert\geq\zeta\right)\leq\frac{C^{\prime\prime}}{\zeta^{2}}\lim_{k\rightarrow\infty}\sum_{i=k}^{\infty}(\frac{\gamma_{i}^{2}}{\delta_{i}^{2}})=0\end{split}.$

The last inequality follows by the assumption
$\sum_{k}(\frac{\gamma_{k}^{2}}{\delta_{k}^{2}})<\infty$. So, by the
convergence of the martingale sequence and Lemma 3, we can conclude that
$x_{k}\to H$ almost surely as $k\to\infty$. ∎

####  4.1.3 Proof of Theorem 2

###### Proof.

This proof follows from Proposition 1 of Chin (1997) after noting the
following facts:

$\displaystyle\mathbb{E}\left[\frac{(u_{k}^{i})^{4}}{(1+\lVert
u\rVert^{2})}\right]$ $\displaystyle\leq\mathbb{E}[\lVert u\rVert^{4}],$
$\displaystyle\mathbb{E}\left[\frac{uu^{T}}{(1+\lVert u\rVert^{2})}\right]$
$\displaystyle\preceq\mathbb{I},$ (22)

∎

####  4.1.4 Proof of Proposition 2

For establishing the avoidance of traps result in Proposition 2, we require a
result from Pemantle (1990), which we state below.

######  Theorem 7.

Let $Y$ be a function in $\Delta\subseteq\mathbb{R}^{d}$ such that
$Y:\Delta\to T\Delta$ where $T\Delta$ is the tangent space of $\Delta$ at each
point. Consider a sequence of random variables $\\{x_{n}:n\geq 0\\}$ that are
updated as in 23 with a given $x_{0}$.

$x_{k+1}=x_{k}+a_{k}Y(x_{k})+\mu_{k}.$ (23)

Let $p$ be any critical point, i.e., $Y(p)=0$, and let $\mathcal{N}$ be a
neighborhood of $p$. Assume that there are constants $\gamma\in(1/2,1]$ and
$c_{5},c_{6},c_{7},c_{8}\geq 0$ for which the following conditions are
satisfied whenever $x_{n}\in\mathcal{N}$ and $n$ is sufficiently large:

  * (i)

$p$ is a linearly unstable critical point of $Y$;

  * (ii)

$\frac{c_{5}}{k^{\gamma}}\leq a_{k}\leq\frac{c_{6}}{k^{\gamma}}$;

  * (iii)

$\mathbb{E}\left[(\mu_{k}\cdot\theta)^{+}|\mathcal{F}_{k}\right]\geq
c_{7}/k^{\gamma}$ for every unit vector $\theta\in T\Delta$;

  * (iv)

$\lVert\mu_{k}\rVert\leq c_{8}/k^{\gamma}$.

In the above,
$(\mu_{k}\cdot\theta)^{+}\stackrel{{\scriptstyle\triangle}}{{=}}\max\\{\mu_{k}\cdot\theta,0\\}$
is the positive part of $\mu_{k}\cdot\theta$. Assume $Y$ is smooth enough (at
least $C^{2}$) to apply the stable manifold theorem. Then $\mathbb{P}(x_{k}\to
p)=0.$

######  Remark 7.

Note that, the limit of $x_{k}$ exists if $Y(x)=0$. Here the iteration rule 23
can be considered as a discrete version of the differential equation
$\dot{x(t)}=Y(x(t))$ with initial condition $x_{0}=v(0)$. We need to classify
the points where $Y(v)=0$. Consider a linear approximation $T$ near a critical
point $p$ of $Y(\cdot)$ such that $Y(p+x)=T(x)+\mathcal{O}(|x|^{2})$. $p$ is
said to be attracting point if real part of the eigenvalue of T is negative
and in such a case $x_{t}$ will converge to $p$ if there are no other
attracting points for the ODE. On the other hand, $p$ is said to be linearly
unstable if some eigenvalue has positive real part and $x_{n}$ exists in the
neighbourhood of $p$ for any choice of $v(0)$ which is not on a stable
manifold of smaller dimension. However, if all the eigenvalues of $T$ have a
positive real part, then $p$ is said to be a repelling node and the sequence
$x_{k},k\geq 0$ will never converge. Theorem 7 gives conditions under which
$\mathbb{P}(x_{k}\to p)=0$ when $p$ is a repelling point as well as a linearly
unstable critical point.

#### Proof of Proposition 2

###### Proof.

Let’s rewrite the update rule 7 as

$x_{k+1}=x_{k}-\gamma_{k}G(x_{k},u_{k},\delta_{k})-\epsilon_{k},$ (24)

where ${\displaystyle
G(x_{k},u_{k},\delta_{k})=\frac{f(x_{k}+\delta_{k}u_{k})-f(x_{k})}{\delta_{k}}\frac{(d+1)u_{k}}{1+\lVert
u_{k}\rVert^{2}}}$ and
${\displaystyle\epsilon_{k}=\gamma_{k}\frac{\eta_{k}^{+}-\eta_{k}}{\delta_{k}}\frac{(d+1)u_{k}}{1+\lVert
u_{k}\rVert^{2}}}$.
Here $\mathbb{E}[\epsilon_{k}|\mathcal{F}_{k}]=0$. Note that 24 is equivalent
to 23 by considering $-\epsilon_{k}=\mu_{k}$ and $-G=Y$.

We will show that the conditions stated in Theorem 7 will hold for our case.
Let $\gamma\in(\frac{1}{2},1]$. Choose the step-size
$\gamma_{k}\in[\frac{c_{8}}{k^{\gamma}},\frac{c_{9}}{k^{\gamma}}]$ such that
both $(ii)$ and (A1) are satisfied, and in (A3), we have considered
$\mathbb{E}|\eta_{k}^{+}|^{2}$ and $\mathbb{E}|\eta_{k}^{+}|^{2}$ are bounded
which in turn implies that $|\eta_{k}^{+}-\eta_{k}|$ is bounded and noting the
fact $u_{k}$ is the truncated Cauchy distribution (over the unit sphere), one
can trivially show that $\lVert\epsilon_{k}\rVert$ is bounded, which implies
$(iv)$.

We now show $(iii)$ for our case below. Note that
$\epsilon_{k}=(\eta_{k}^{+}-\eta_{k})u_{k}m_{k}$, where
$m_{k}=\frac{(d+1)\gamma_{k}}{1+\lVert u_{k}\rVert^{2}}$. Consider the unit
vector with the $i$th entry as $1$, i.e., $\theta=(0,0,...1,..0)^{T}$. This
implies $\epsilon_{k}\cdot\theta=(\eta_{k}^{+}-\eta_{k})u_{k}^{i}m_{k}$. Here,
$u_{k}^{i}$ denotes the $i$th entry of the vector $u$ at the $k$th iteration.

$\displaystyle\mathbb{E}[(\epsilon_{k}\cdot\theta)^{+}]$
$\displaystyle=\mathbb{E}[((\eta_{k}^{+}-\eta_{k})\cdot u_{k}^{i}m_{k})^{+}]$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{=}}\mathbb{E}[\frac{(\eta_{k}^{+}-\eta_{k})\cdot
u_{k}^{i}m_{k}+|(\eta_{k}^{+}-\eta_{k})\cdot u_{k}^{i}m_{k}|}{2}]$
$\displaystyle\stackrel{{\scriptstyle(c)}}{{=}}\mathbb{E}[\frac{|(\eta_{k}^{+}-\eta_{k})\cdot
u_{k}^{i}m_{k}|}{2}]$
$\displaystyle=\mathbb{E}[\frac{|\eta_{k}^{+}-\eta_{k}|\cdot|u_{k}^{i}m_{k}|}{2}]$
$\displaystyle\stackrel{{\scriptstyle(d)}}{{\geq}}\frac{c_{9}c_{10}}{2}.$

In the above, we used (i) the fact that $\max(x,y)=\frac{x+y+|x-y|}{2}$ to
infer the equality in $(b)$; (ii) $\mathbb{E}[(\eta_{k}^{+}-\eta_{k})\cdot
u_{k}^{i}m_{k}]=0$ since $E[\eta_{k}^{+}-\eta_{k}\mid\mathcal{F}_{k}]=0$ and
$u_{k}$ is independent of $\mathcal{F}_{k}$, to infer the equality in $(c)$;
and (iii) (A6) in conjunction with $\mathbb{E}_{u}|u_{k}^{i}m_{k}|\geq c_{10}$
for some positive constant $c_{10}$ is used to infer $(c)$. Thus, condition
(iii) of Theorem 7 holds. The main claim now follows by an application of
Theorem 7.

∎

###  4.2 Non-asymptotic convergence proofs

####  4.2.1 Proof of Theorem 3

###### Proof.

We use the the proof technique from Bhavsar and Prashanth (2022)(in
particular, proposition-1 there) in order to proof the main claim here.
However unlike Bhavsar and Prashanth (2022) we have a gradient estimate that
comes from truncated Cauchy distribution.
Lets define $\alpha_{k}\equiv(\xi_{k},\xi_{k}^{+},u_{k},\delta),k\geq 1$ and
$\alpha_{[N]}:=(\alpha_{1},\alpha_{2},...\alpha_{N})$ U̇sing Taylor series
expansion over $f(x_{k})$ for any $k=1,2,...,N$ the following is obtained

$\begin{split}f(x_{k+1})&\leq f(x_{k})-\gamma_{k}\langle\nabla
f(x_{k}),G(x_{k},\alpha)\rangle+\frac{L}{2}\gamma_{k}^{2}\left\lVert
G(x_{k},\alpha_{k})\right\rVert^{2},\\\
&=f(x_{k})-c_{2}\gamma_{k}\left\lVert\nabla
f(x_{k})\right\rVert^{2}-\gamma_{k}\langle\nabla
f(x_{k}),\Gamma_{k}\rangle+\frac{L}{2}\gamma_{k}^{2}\left\lVert
G(x_{k},\alpha_{k})\right\rVert^{2}.\end{split}$

Here $\Gamma_{k}\equiv G(x_{k},\xi_{k},u_{k},\delta)-c_{2}\nabla f(x)\equiv
G(x_{k},\alpha_{k})-c_{2}\nabla f(x_{k})$. Adding upto $N$-terms both side of
these inequalities and considering $f^{*}\leq f(x_{N+1})$ and
$\gamma_{k}=\gamma$ for all $k$, we obtain

$\displaystyle\sum_{k=1}^{N}c_{2}\gamma\left\lVert\nabla
f(x_{k})\right\rVert^{2}$ $\displaystyle\leq
f(x_{1})-f^{*}-\sum_{k=1}^{N}\gamma\langle\nabla
f(x_{k}),\gamma\rangle+\frac{L}{2}\sum_{k=1}^{N}\gamma^{2}\left\lVert
G(x_{k},\alpha_{k})\right\rVert^{2}.$ (25)

Now by Lemma 1 we have

$\mathbb{E}_{\alpha_{[k]}}[\Gamma_{k}]=\mathbb{E}_{\alpha_{[k]}}[\Gamma_{k}|\alpha_{k-1}]=\mathbb{E}_{\alpha_{[k]}}[\Gamma_{k}|x_{k}]=\mathbb{E}_{\alpha[k]}[G-c_{2}\nabla
f|x_{k}]\stackrel{{\scriptstyle(e)}}{{\leq}}\tau\mathbf{1}_{d\times 1},$

notice $\tau=c_{2}^{\prime}\delta$ where $c_{2}^{\prime}=(d+1)B$ is a constant
that arise from the Taylor series as defined. In the above vector inequality
$(e)$ is element-wise. From Section 4.1 we have

$\mathbb{E}_{\alpha_{[k]}}[\lVert
G\rVert^{2}]\leq\lVert\mathbb{E}_{\alpha_{[k]}}[G]\rVert^{2}+\frac{C^{\prime\prime}}{\delta^{2}}.$

Hence by taking the expectation w.r.t $\alpha_{[N]}$ on both side of 25 the
following is obtained

$\displaystyle\sum_{k=1}^{N}c_{2}\gamma\mathbb{E}_{\alpha_{[N]}}\left\lVert\nabla
f(x_{k})\right\rVert^{2}$ $\displaystyle\leq
D+BNc_{2}^{\prime}\delta\gamma+\frac{L}{2}\sum_{k=1}^{N}\gamma^{2}\Big{[}\mathbb{E}_{\alpha_{[N]}}\lVert\nabla
f(x_{k})\lVert^{2}+2c_{2}^{\prime}\delta B+dc_{2}^{\prime
2}\delta^{2}+\frac{C^{\prime\prime}}{\delta^{2}}\Big{]}.$

The above inequality uses the fact $-\lVert
V\rVert_{1}\leq\sum_{k=1}^{d}v_{k}$ for a $d-$dimensional vector $V$ followed
by $\lVert\nabla f(x_{k})\rVert_{1}\leq\lVert\nabla f(x_{k})\rVert\leq B$ from
(A7). Note that $D=f(x_{1})-f^{*}$. By rearranging the terms we have

$\displaystyle\left[c_{2}\gamma-\frac{L\gamma^{2}}{2}\right]\sum_{k=1}^{N}\mathbb{E}_{\alpha_{[N]}}\left\lVert\nabla
f(x_{k})\right\rVert^{2}$ $\displaystyle\leq
D+BNc_{2}^{\prime}\delta\gamma+\frac{LN}{2}\left(2c_{2}^{\prime}\delta
B+dc_{2}^{\prime
2}\delta^{2}+\frac{C^{\prime\prime}}{\delta^{2}}\right)\gamma^{2}.$

Due to the choice of
$\gamma_{k}=\gamma=\min\left\\{\frac{c_{2}}{L},\frac{1}{N^{2/3}}\right\\}$, it
is obvious that $N\left[c_{2}\gamma-\frac{L\gamma^{2}}{2}\right]\geq 0$. Thus
by dividing both sides of the above inequality by
$N\left[c_{2}\gamma-\frac{L\gamma^{2}}{2}\right]$ and noting the fact

$\mathbb{P}_{R}(k)=Prob(R=k)=\frac{\left[c_{2}\gamma-\frac{L\gamma^{2}}{2}\right]}{N\left[c_{2}\gamma-\frac{L\gamma^{2}}{2}\right]}=\frac{1}{N},$
(26)

the following is obtained

$\displaystyle\mathbb{E}_{\alpha_{[N]}}\left\lVert\nabla
f(x_{R})\right\rVert^{2}$
$\displaystyle\leq\frac{1}{N\left[c_{2}\gamma-\frac{L\gamma^{2}}{2}\right]}\bigg{[}D+BNc_{2}^{\prime}\delta\gamma+\frac{LN}{2}\left(2c_{2}^{\prime}\delta
B+dc_{2}^{\prime
2}\delta^{2}+\frac{C^{\prime\prime}}{\delta^{2}}\right)\gamma^{2}\bigg{]}.$

By considering, $\gamma=\frac{c_{2}}{L}$, the following is obtained

$Nc_{2}\gamma\left[1-\frac{L}{2c_{2}}\gamma\right]\geq\frac{Nc_{2}\gamma}{2}.$

From the above inequality, we can write

$\displaystyle\mathbb{E}\left[\left\lVert\nabla
f(x_{R})\right\rVert^{2}\right]$
$\displaystyle\leq\frac{2D}{Nc_{2}\gamma}+\frac{2Bc_{2}^{\prime}\delta}{c_{2}}+\frac{L}{c_{2}}\left(2c_{2}^{\prime}\delta
B+dc_{2}^{\prime
2}\delta^{2}+\frac{C^{\prime\prime}}{\delta^{2}}\right)\gamma$
$\displaystyle\leq\frac{2D}{Nc_{2}}\max\left\\{\frac{L}{c_{2}},N^{2/3}\right\\}+\frac{2Bc_{2}^{\prime}\delta}{c_{2}}+\frac{L}{c_{2}N^{2/3}}\left(2c_{2}^{\prime}\delta
B+dc_{2}^{\prime 2}\delta^{2}+\frac{C^{\prime\prime}}{\delta^{2}}\right)$
$\displaystyle\stackrel{{\scriptstyle(g)}}{{\leq}}\Big{(}\frac{2DL}{Nc_{2}^{2}}+\frac{2D}{N^{1/3}}\Big{)}+\frac{2Bc_{2}^{\prime}}{c_{2}N^{1/6}}+\frac{L}{c_{2}N^{2/3}}\left(\frac{2Bc_{2}^{\prime}}{N^{1/6}}+\frac{dc_{2}^{\prime
2}}{N^{1/3}}+\frac{C^{\prime\prime}}{N^{-1/3}}\right)$
$\displaystyle=\Big{(}\frac{2DL}{c_{2}N}+\frac{2D}{N^{1/3}}\Big{)}+\frac{2Bc_{2}^{\prime}}{c_{2}N^{1/6}}+\frac{L}{c_{2}}\left(\frac{2Bc_{2}^{\prime}}{N^{5/6}}+\frac{dc_{2}^{\prime
2}}{N}+\frac{C^{\prime\prime}}{N^{1/3}}\right).$

Note that $(g)$ uses the condition of $\delta=\frac{1}{N^{1/6}}$. Hence
proved. ∎

####  4.2.2 Proof of Theorem 5

The proof is accumulated from a sequence of lemmas..We follow the technique
from Ghadimi and Lan (2013), and our proof of the lemmas involves significant
deviations owing to the fact that biased gradient information under truncated
Cauchy distribution are available instead of unbiased gradient information
(under Gaussian distribution).

######  Lemma 4.

let $f$ be the function satisfying (A7) and (A8) and $\delta$ be the smoothing
parameter defined in Algorithm 1. Then

$\lVert
f_{\delta}(x)\rVert\leq\frac{c_{11}}{d+1}\left[\frac{L\delta^{2}}{2}+2\delta
B\right].$

###### Proof.

$\begin{split}f_{\delta}(x)&=\mathbb{E}_{u}(f(x+\delta u)-f(x))\\\
&=\mathbb{E}_{u}[f(x+\delta u)-f(x)-\delta\langle\nabla
f(x),u\rangle]+\delta\mathbb{E}_{u}[\langle\nabla f(x),u\rangle]\\\
|f_{\delta}(x)|&\stackrel{{\scriptstyle(d)}}{{\leq}}\frac{L\delta^{2}}{2}\mathbb{E}_{u}[\lVert
u\rVert^{2}]+\delta\lVert\nabla f(x)\rVert\mathbb{E}_{u}[\lVert u\rVert]\\\
&\stackrel{{\scriptstyle(d)}}{{\leq}}\frac{\delta^{2}c_{11}L}{2(d+1)}+\frac{2\delta
Bc_{11}}{2d+1}\\\ &\leq\frac{c_{11}}{d+1}\left[\frac{L\delta^{2}}{2}+2\delta
B\right].\end{split}$

where $(d)$ follows from smoothness of $f$ and Cauchy-Schwarz inequality and
$(e)$ follows from Lemma 2. ∎

######  Proposition 3.

The solution of a linear system of equations $Ax=y$, where $A$ is a non-
invertible matrix, is $Py$ where $P$ is the generalized inverse of the matrix
$A$ that satisfies $APA=A$ and $PAP=P$.

######  Lemma 5.

Under (A7) and (A8), we have

  1. (a)

$\lVert\nabla f(x)\rVert^{2}\leq 2c_{12}\lVert\nabla
f_{\delta}(x)\rVert^{2}+\frac{c_{11}c_{12}\delta^{2}L^{2}(d+1)}{2}.$

  2. (b)

$\lVert\nabla
f_{\delta}(x)\rVert^{2}\leq\frac{2c_{11}}{d+1}\left[\frac{(d+1)^{2}\delta^{2}L^{2}}{4}+\lVert\nabla
f(x)\rVert^{2}\right].$

###### Proof.

Notice that

$\begin{split}\mathbb{E}_{u}\left[\frac{(d+1)u}{1+\lVert
u\rVert^{2}}\langle\nabla f(x),u\rangle\right]&=A\nabla f(x)\\\ \end{split}$
(27)

where $A=\mathbb{E}_{u}\left(\frac{(d+1)uu^{T}}{1+\lVert u\rVert^{2}}\right)$
is a matrix. Let $P$ be the generalized inverse of $A$ which satisfies the
condition in Proposition 3. Now

$\displaystyle\nabla f(x)$
$\displaystyle\stackrel{{\scriptstyle(e)}}{{=}}P\cdot\mathbb{E}\bigg{[}\frac{(d+1)u}{1+\lVert
u\rVert^{2}}\langle\nabla f(x),u\rangle\bigg{]},$
$\displaystyle=\frac{P}{\delta}\mathbb{E}\left[\frac{(d+1)u}{1+\lVert
u\rVert^{2}}\\{f(x+\delta
u)-f(x)\\}\right]-\frac{P}{\delta}\mathbb{E}\bigg{[}\frac{(d+1)u}{1+\lVert
u\rVert^{2}}\\{f(x+\delta u)-f(x)-\delta\langle\nabla
f(x),u\rangle\\}\bigg{]},$

$(e)$ is a simple application of Proposition 3. Let $\lVert\cdot\rVert_{F}$
denotes the Frobenius norm of the matrix. Then $\lVert Px\rVert_{2}\leq
c_{12}\lVert x\rVert_{2}$ where $c_{12}=\lVert P\rVert_{F}$

$\displaystyle\lVert\nabla f(x)\rVert^{2}$
$\displaystyle\stackrel{{\scriptstyle(f)}}{{\leq}}2\lVert
P\rVert_{F}^{2}\lVert\nabla f_{\delta}(x)\rVert^{2}+\frac{2\lVert
P\rVert_{F}^{2}}{\delta^{2}}\mathbb{E}\left[\frac{(d+1)^{2}\lVert
u\rVert^{6}L^{2}\delta^{4}}{4(1+\lVert u\rVert^{2})^{2}}\right],$
$\displaystyle=2c_{12}\lVert\nabla
f_{\delta}(x)\rVert^{2}+\frac{c_{12}\delta^{2}L^{2}(d+1)^{2}}{2}\mathbb{E}\left[\frac{\lVert
u\rVert^{6}}{(1+\lVert u\rVert^{2})^{2}}\right],$ $\displaystyle\leq
2c_{12}\lVert\nabla
f_{\delta}(x)\rVert^{2}+\frac{c_{12}\delta^{2}L^{2}(d+1)^{2}}{2}\mathbb{E}\left[\lVert
u\rVert^{6}\right],$ $\displaystyle\leq 2c_{12}\lVert\nabla
f_{\delta}(x)\rVert^{2}+\frac{c_{12}\delta^{2}L^{2}(d+1)^{2}}{2}\mathbb{E}\left[\lVert
u\rVert^{6}\right]$ $\displaystyle\leq 2c_{12}\lVert\nabla
f_{\delta}(x)\rVert^{2}+\frac{c_{11}c_{12}\delta^{2}L^{2}(d+1)^{2}}{2(d+3)},$
$\displaystyle\leq 2c_{12}\lVert\nabla
f_{\delta}(x)\rVert^{2}+c_{11}c_{12}\delta^{2}L^{2}(d+1),$

where $(f)$ uses the Cauchy-Schwarz inequality followed by Lemma 2.

We now prove the second claim.

$\begin{split}\left[f(x+\delta u)-f(x)\right]^{2}&=[f(x+\delta
u)-f(x)-\delta\langle\nabla f(x),u\rangle+\delta\langle\nabla
f(x),u\rangle]^{2},\\\ &\leq 2\left(\frac{\delta^{2}}{2}L\left\lVert
u\right\lVert^{2}\right)^{2}+2\delta^{2}\langle\nabla f(x),u\rangle^{2},\\\
&\leq 2\left(\frac{\delta^{2}}{2}L\left\lVert
u\right\lVert^{2}\right)^{2}+2\delta^{2}\lVert\nabla f(x)\rVert^{2}\left\lVert
u\right\rVert^{2}.\\\ \end{split}$

Now, by taking the norm in the both side of the inequality and employing Lemma
2 the following is obtained

$\displaystyle\lVert\nabla
f_{\delta}(x)\rVert^{2}\leq\frac{1}{\delta^{2}}\mathbb{E}_{u}\left[\left(f(x+\delta
u)-f(x)\right)^{2}\frac{(d+1)^{2}\left\lVert u\right\lVert^{2}}{(1+\lVert
u\rVert^{2})^{2}}\right],$ $\displaystyle\leq\frac{((d+1)\delta
L)^{2}}{2}\mathbb{E}_{u}\left[\frac{\left\lVert u\right\lVert^{6}}{(1+\lVert
u\rVert^{2})^{2}}\right]+2\lVert\nabla
f(x)\rVert^{2}(d+1)^{2}\mathbb{E}_{u}\left[\frac{\left\lVert
u\right\lVert^{4}}{(1+\lVert u\rVert^{2})^{2}}\right],$
$\displaystyle\leq\frac{((d+1)\delta L)^{2}}{2}\mathbb{E}_{u}\left\lVert
u\right\lVert^{6}+2\lVert\nabla
f(x)\rVert^{2}(d+1)^{2}\mathbb{E}_{u}\left\lVert u\right\lVert^{4},$
$\displaystyle\leq\frac{((d+1)\delta L)^{2}c_{11}}{2(d+3)}+\frac{2\lVert\nabla
f(x)\rVert^{2}c_{11}}{d+2},$
$\displaystyle\leq\frac{2c_{11}}{d+1}\left[\frac{(d+1)^{2}\delta^{2}L^{2}}{4}+\lVert\nabla
f(x)\rVert^{2}\right].$

∎

In the following lemma we state and prove a general result that holds true for
any choice of non-increasing stepsize sequence, smoothing parameter.
Subsequently, we specialize the result for the choice of parameters suggested
in Theorem 5, to prove the same.

######  Lemma 6.

Suppose, $x_{k}$ generated by Algorithm 1 and $\\{\gamma_{k}\\}$ be the
desired step-sizes for the iteration. Let and the probability mass function
$P_{R}(\cdot)$ are chosen such that $\gamma_{k}\leq\frac{1}{2Lc_{13}}$ and

$P_{R}(k):=Prob(R=k)=\frac{\left[\gamma_{k}-Lc_{13}\gamma_{k}^{2}\right]}{\sum_{k=1}^{N}\left[\gamma_{k}-Lc_{13}\gamma_{k}^{2}\right]},$
(28)

where $c_{13}=\frac{4c_{11}c_{12}}{d+1}.$ and $c_{11}$,$c_{12}$ are defined in
16 ,27 respectively. Then under (A7),(A9) and for any $N\geq 1$

$\displaystyle\mathbb{E}\left[\left\lVert\nabla
f(x_{R})\right\rVert^{2}\right]$
$\displaystyle\leq\frac{1}{\sum_{k=1}^{N}\left[\gamma_{k}-Lc_{13}\gamma_{k}^{2}\right]}\bigg{[}\frac{L\delta^{2}c_{13}}{2}+2c_{13}\delta
B+L^{2}\delta^{2}d^{2}c_{13}\sum_{k}\gamma_{k}$ (29)
$\displaystyle+Lc_{13}\left(\frac{d\delta^{2}L^{2}}{4}+\sigma^{2}\right)\sum\gamma_{k}^{2}\bigg{]}.$
(30)

###### Proof.

Let $\alpha_{k}\equiv(\xi_{k},u_{k}),k\geq 1$ and
$\alpha_{[N]}:=(\alpha_{1},\alpha_{2},...\alpha_{N})$. Denote
$\Gamma_{k}\equiv G_{\delta}(x_{k},\xi_{k},u_{k})-\nabla f_{\delta}(x)\equiv
G_{\delta}(x_{k},\alpha_{k})-\nabla f_{\delta}(x)$. Now, by Taylor series
expansion $f_{\delta}(x_{k+1})$ we have for any $k=1,2,...,N$,

$\begin{split}f_{\delta}(x_{k+1})&\leq
f_{\delta}(x_{k})-\gamma_{k}\langle\nabla
f_{\delta}(x_{k}),G_{\delta}(x_{k},\alpha)\rangle+\frac{L}{2}\gamma_{k}^{2}\left\lVert
G_{\delta}(x_{k},\alpha_{k})\right\rVert^{2}\\\
&=f_{\delta}(x_{k})-\gamma_{k}\left\lVert\nabla
f_{\delta}(x_{k})\right\rVert^{2}-\gamma_{k}\langle\nabla
f_{\delta}(x_{k}),\Gamma_{k}\rangle+\frac{L}{2}\gamma_{k}^{2}\left\lVert
G_{\delta}(x_{k},\alpha_{k})\right\rVert^{2}.\end{split}$

Adding upto $N$-terms both side of these inequalities and applying
$f_{\delta}^{*}\leq f_{\delta}(x_{N+1})$, we obtain

$\displaystyle\sum_{k=1}^{N}\gamma_{k}\left\lVert\nabla
f_{\delta}(x_{k})\right\rVert^{2}$ $\displaystyle\leq
f_{\delta}(x_{1})-f_{\delta}^{*}-\sum_{k=1}^{N}\gamma_{k}\langle\nabla
f_{\delta}(x_{k}),\Gamma_{k}\rangle+\frac{L}{2}\sum_{k=1}^{N}\gamma_{k}^{2}\left\lVert
G_{\delta}(x_{k},\alpha_{k})\right\rVert^{2}.$ (31)

From the unbiased property of $G_{\delta}(x_{k},u_{k},\xi_{k})$ the following
holds

$\mathbb{E}[\langle\nabla
f_{\delta}(x_{k}),\Gamma_{k}\rangle|\alpha_{[k-1]}]=0.$

Now by (A9) and Lemma 5(b)

$\begin{split}\mathbb{E}[\left\lVert
G_{\delta}(x_{k},\alpha_{k})\right\rVert^{2}|\alpha_{[k-1]}]&\leq\frac{2c_{11}}{d+1}\left[\mathbb{E}\lVert\nabla
F(x_{k},\xi_{k})\rVert^{2}|\alpha_{[k-1]}]+\frac{(d+1)^{2}\delta^{2}L^{2}}{4}\right]\\\
&\leq\frac{2c_{11}}{d+1}\left[2\mathbb{E}[\left\lVert\nabla
f(x_{k})\right\rVert^{2}|\alpha_{[k-1]}+\sigma^{2}]+\frac{(d+1)^{2}\delta^{2}L^{2}}{4}\right].\end{split}$
(32)

Note that,the second inequality implies from variance bound of (A7). Taking
expectations with respect to $\alpha_{[N]}$ on both sides of 31, we have

$\displaystyle\sum_{k=1}^{N}\gamma_{k}\mathbb{E}_{\alpha_{[N]}}\left\lVert\nabla
f_{\delta}(x_{k})\right\rVert^{2}$ $\displaystyle\leq
f_{\delta}(x_{1})-f_{\delta}^{*}+\frac{L}{2}\sum_{k=1}^{N}\gamma_{k}^{2}\bigg{[}\frac{2c_{11}}{d+1}\Big{[}2\mathbb{E}[\left\lVert\nabla
f(x_{k})\right\rVert^{2}|\alpha_{[k-1]}+\sigma^{2}]+\frac{(d+1)^{2}\delta^{2}L^{2}}{4}\Big{]}\bigg{]}.$

Applying Lemma 5(a) in the r.h.s and rearranging the term we obtain

$\displaystyle\sum_{k=1}^{N}\gamma_{k}\left[\frac{1}{2c_{12}}\mathbb{E}_{\alpha_{[N]}}\left\lVert\nabla
f(x_{k})\right\rVert^{2}-\frac{c_{11}(d+1)L^{2}\delta^{2}}{2}\right]$
$\displaystyle\leq
f_{\delta}(x_{1})-f_{\delta}^{*}+\frac{LNc_{11}}{d+1}\left(\frac{(d+1)\delta^{2}L^{2}}{4}+2\sigma^{2}\right)\sum_{k=1}^{N}\gamma_{k}^{2}+\frac{2Lc_{11}}{d+1}\sum_{k=1}^{N}\gamma_{k}^{2}\mathbb{E}[\lVert\nabla
f(x_{k})\rVert^{2}].$

From Lemma 4 and using the fact
$f_{\delta}(x_{1})-f_{\delta}^{*}\leq|f_{\delta}(x_{1})-f_{\delta}^{*}|\leq|f_{\delta}(x_{1})|+|f_{\delta}^{*}|$
one can have

$\displaystyle\sum_{k=1}^{N}\left[\frac{\gamma_{k}}{2c_{12}}-\frac{2Lc_{11}\gamma_{k}^{2}}{d+1}\right]\mathbb{E}_{\alpha_{[N]}}\left\lVert\nabla
f(x_{k})\right\rVert^{2}$
$\displaystyle\leq\frac{c_{11}L\delta^{2}}{d+1}+\frac{4c_{11}\delta
B}{d+1}+\frac{NL^{2}\delta^{2}(d+1)c_{11}}{2}\sum_{k=1}^{N}\gamma_{k}+\frac{LNc_{11}}{d+1}\left(\frac{(d+1)\delta^{2}L^{2}}{4}+2\sigma^{2}\right)\sum_{k=1}^{N}\gamma_{k}^{2}.$

Thus by rearranging the above inequality we obtain

$\displaystyle\sum_{k=1}^{N}\left[\gamma_{k}-\frac{4c_{11}c_{12}L\gamma_{k}^{2}}{d+1}\right]\mathbb{E}_{\alpha_{[N]}}\left\lVert\nabla
f(x_{k})\right\rVert^{2}$
$\displaystyle\leq\frac{2c_{11}c_{12}L\delta^{2}}{d+1}+\frac{8c_{11}c_{12}\delta
B}{d+1}+L^{2}\delta^{2}(d+1)c_{11}c_{12}\sum_{k=1}^{N}\gamma_{k}+\frac{2Lc_{11}c_{12}}{d+1}\left(\frac{(d+1)\delta^{2}L^{2}}{4}+2\sigma^{2}\right)\sum_{k=1}^{N}\gamma_{k}^{2}.$

Let $c_{13}=\frac{4c_{11}c_{12}}{d+1}$

$\displaystyle\sum_{k=1}^{N}\left[\gamma_{k}-Lc_{13}\gamma_{k}^{2}\right]\mathbb{E}_{\alpha_{[N]}}\left\lVert\nabla
f(x_{k})\right\rVert^{2}$ (33)
$\displaystyle\leq\frac{L\delta^{2}c_{13}}{2}+2c_{13}\delta
B+\frac{L^{2}\delta^{2}(d+1)^{2}c_{13}}{4}\sum_{k=1}^{N}\gamma_{k}+\frac{Lc_{13}}{2}\left(\frac{(d+1)\delta^{2}L^{2}}{4}+\sigma^{2}\right)\sum_{k=1}^{N}\gamma_{k}^{2}.$
(34)

Dividing both sides of the above inequality by
$\sum_{k=1}^{N}\left[\gamma_{k}-Lc_{13}\gamma_{k}^{2}\right]$ and notice that

$\begin{split}\mathbb{E}\left[\left\lVert\nabla
f(x_{R})\right\rVert^{2}\right]&=\mathbb{E}_{R,\alpha_{[N]}}\left[\left\lVert\nabla
f(x_{R})\right\rVert^{2}\right]=\frac{\sum_{k=1}^{N}\left[\gamma_{k}-Lc_{13}\gamma_{k}^{2}\right]\mathbb{E}_{\alpha_{[N]}}\left\lVert\nabla
f(x_{k})\right\rVert^{2}}{\sum_{k=1}^{N}\left[\gamma_{k}-Lc_{13}\gamma_{k}^{2}\right]},\end{split}$
(35)

we obtain 29 by replacing $d+1$ with $2d$ in 33.

∎

We now specialize the result obtained from Lemma 6 to get the tight bound for
Algorithm 1 as describe in Theorem 5.

#### Proof of Theorem 5

###### Proof.

Recall the step-size $\gamma_{k}=\gamma$, smoothing parameter
$\delta_{k}=\delta,\forall k\geq 1,$ where

$\begin{split}\gamma_{k}=\gamma&=\min\left\\{\frac{1}{2Lc_{13}},\frac{1}{c_{13}\sigma\sqrt{N}}\right\\}\
,\delta=\frac{1}{L\sqrt{dNc_{13}}}.\end{split}$ (36)

From the above condition of $\gamma$ by considering
$\gamma\leq\frac{1}{2Lc_{13}}$ we have

$N\gamma\left[1-Lc_{13}\gamma\right]\geq\frac{N\gamma}{2}.$ (37)

Henceforth, from the above inequality and by Lemma 6, we obtain

$\displaystyle\mathbb{E}\left[\left\lVert\nabla
f(x_{R})\right\rVert^{2}\right]$
$\displaystyle\stackrel{{\scriptstyle(h)}}{{\leq}}\frac{L\delta^{2}c_{13}}{N\gamma}+\frac{4\delta
Bc_{13}}{N\gamma}+2L^{2}\delta^{2}d^{2}c_{13}+2Lc_{13}\left(\frac{d\delta^{2}L^{2}}{4}+\sigma^{2}\right)\gamma,$
$\displaystyle\leq\left(\frac{L\delta^{2}c_{13}}{N}+\frac{4\delta
Bc_{13}}{N}\right)\max\\{2Lc_{13},\sigma
c_{13}\sqrt{N}\\}+2L^{2}\delta^{2}d^{2}c_{13}+\frac{d\delta^{2}L^{2}}{4}+\frac{2\sigma
L}{\sqrt{N}},$
$\displaystyle\leq\left(\frac{L\delta^{2}c_{13}}{N}+\frac{4\delta
Bc_{13}}{N}\right)(2Lc_{13}+\sigma
c_{13}\sqrt{N})+L^{2}\delta^{2}d(2c_{13}d+1)+\frac{4\sigma L}{\sqrt{N}},$
$\displaystyle\stackrel{{\scriptstyle(k)}}{{=}}\frac{c_{13}^{2}}{N}(L\delta^{2}+4\delta
B)(2L+\sigma\sqrt{N})+L^{2}\delta^{2}d(2c_{13}d+1)+\frac{2\sigma
L}{\sqrt{N}},$
$\displaystyle\leq\frac{c_{13}}{LN}(\frac{1}{Nd}+\frac{4B\sqrt{c_{13}}}{\sqrt{Nd}})(2L+\sigma\sqrt{N})+\frac{(2c_{13}d+1)}{Nc_{13}}+\frac{2\sigma
L}{\sqrt{N}},$
$\displaystyle=\frac{c_{13}}{LN}(\frac{2L}{Nd}+\frac{\sigma}{d\sqrt{N}}+\frac{8LB\sqrt{c_{13}}}{\sqrt{dN}}+\frac{\sigma
B\sqrt{c_{13}}}{\sqrt{d}})+\frac{(2c_{13}d+1)}{Nc_{13}}+\frac{2\sigma
L}{\sqrt{N}},$
$\displaystyle\leq\frac{c_{13}}{LN}(2L+\sigma+8LB\sqrt{c_{13}}+\sigma
B\sqrt{c_{13}})+\frac{(2c_{13}d+1)}{Nc_{13}}+\frac{2\sigma L}{\sqrt{N}},$
$\displaystyle=\frac{c_{14}}{N}+\frac{4\sigma L}{\sqrt{N}}.$

Where $c_{14}=2L+c_{13}\sigma+8LBc_{13}^{3/2}+\frac{\sigma
Bc_{13}^{3/2}B}{L}+\frac{2c_{13}d+1}{L}$. Notice that, $(h)$ uses the
condition of $\gamma_{k}$ in 36 and $(k)$ follows from the condition of
$\delta=\frac{1}{L\sqrt{dNc_{13}}}$. Thus by rearranging the terms as noting
the fact that $1/Nd\leq 1$ we get the convergence rate
$\mathcal{O}(1/\sqrt{N}).$ ∎

###  4.3 Non-asymptotic analysis for the balanced estimator

Till now we have covered the analysis with imbalance gradient estimator.
However, with this estimator we will get high bias which is
$\mathcal{O}(\delta_{k})$. We will now introduce a balanced gradient estimator
and later we will show how it will help to achieve low bias as compared to the
previous one.
Recall the finite difference gradient estimate is defined as

$\begin{split}\nabla
f_{\delta}(x)=\frac{1}{\delta}\mathbb{E}_{h(u)}\left[(f(x+\delta
u)-f(x))\frac{(d+1)u}{(1+\lVert u\rVert^{2})}\right].\\\ \end{split}$ (38)

Note that by change of variable we can rewrite 38 as

$\begin{split}\nabla
f_{\delta}(x)=\frac{1}{\delta}\mathbb{E}_{h(u)}\left[(f(x)-f(x-\delta
u))\frac{(d+1)u}{(1+\lVert u\rVert^{2})}\right].\\\ \end{split}$ (39)

Now by summing up 38 and 39 one can get the balanced estimator as

$\begin{split}\nabla
f_{\delta}(x)=\frac{1}{2\delta}\mathbb{E}_{h(u)}\left[(f(x+\delta
u)-f(x-\delta u))\frac{(d+1)u}{(1+\lVert u\rVert^{2})}\right].\\\ \end{split}$
(40)

Thus for the case of noisy function measurements the balanced estimator is

$\displaystyle\tilde{G}(x_{k},\xi_{k}^{+},\xi_{k}^{-},u_{k},\delta_{k})$
$\displaystyle\stackrel{{\scriptstyle\triangle}}{{=}}{\displaystyle\bigg{(}\frac{F(x_{k}+\delta_{k}u_{k},\xi_{k}^{+})-F(x_{k}-\delta_{k}u_{k},\xi_{k}^{-})}{2\delta_{k}}\bigg{)}\frac{(d+1)u_{k}}{1+\lVert
u_{k}\rVert^{2}}}.$ (41)

####  4.3.1 Proof of Theorem 4

######  Lemma 7.

Under (A1)-(A4), and $\tilde{G}$ defined in 13 we have almost surely

$\displaystyle\mathbb{E}[\tilde{G}(x_{k},\xi_{k},\xi_{k}^{+},u_{k},\delta_{k})|\mathcal{F}_{k}]=c_{2}\nabla
f(x_{k})+c_{2}^{\prime\prime}\delta^{2}_{k}\mathbf{1}_{d},$ (42)

where $c_{2}=\mathbb{E}_{h(u)}\left[\frac{(d+1)(u_{k}^{i})^{2}}{1+\lVert
u_{k}\rVert^{2}}\right]$, $c_{2}^{\prime\prime}=\frac{B_{1}d^{4}}{3}$ and
$\mathbf{1}_{d}$ is the $d$-dimensional vector of all ones.

###### Proof.

Using Taylor series expansion for truncated Cauchy perturbations we obtain:
$f(x_{k}+\delta_{k}u_{k})-f(x_{k}-\delta_{k}u_{k})$ $=$
$2\delta_{k}u_{k}^{T}\nabla
f(x_{k})+\frac{\delta_{k}^{3}}{6}(\nabla^{3}f(\bar{x_{k}}^{+})+\nabla^{3}f(\bar{x_{k}}^{-}))(u_{k}\otimes
u_{k}\otimes u_{k})$.
Here $\otimes$ denotes the Kronecker product and $\bar{x}^{+}$(respectively,
$\bar{x}^{-}$) are on the line segment between $x$ and $x+\delta u$
(respectively, $x-\delta u)$. So

$\displaystyle\mathbb{E}[\tilde{G}(x_{k},\xi_{k},\xi_{k}^{+},u_{k},\delta_{k})|\mathcal{F}_{k}]$
$\displaystyle=\mathbb{E}\Big{[}\frac{F(x+\delta_{k}u_{k},\epsilon_{k}^{+})-F(x-\delta_{k}u_{k},\epsilon_{k}^{-})}{2\delta_{k}}\frac{(d+1)u_{k}}{1+\lVert
u_{k}\rVert^{2}}|\mathcal{F}_{k}\Big{]}$
$\displaystyle=\mathbb{E}\Big{[}\frac{f(x+\delta_{k}u_{k})-f(x-\delta_{k}u_{k})}{2\delta_{k}}\frac{(d+1)u_{k}}{1+\lVert
u_{k}\rVert^{2}}|\mathcal{F}_{k}\Big{]}+\mathbb{E}\Big{[}\frac{\eta_{k}^{+}-\eta_{k}^{-}}{2\delta_{k}}\frac{(d+1)u_{k}}{1+\lVert
u_{k}\rVert^{2}}|\mathcal{F}_{k}\Big{]},$
$\displaystyle=\mathbb{E}\Big{[}\frac{(d+1)u_{k}u_{k}^{T}\nabla
f(x_{k})}{1+\lVert
u_{k}\rVert^{2}}|\mathcal{F}_{k}\Big{]}+\mathbb{E}\Big{[}\frac{\delta_{k}^{2}(d+1)u_{k}}{12(1+\lVert
u_{k}\rVert^{2})}(\nabla^{3}f(\bar{x_{k}}^{+})+\nabla^{3}f(\bar{x_{k}}^{-}))(u_{k}\otimes
u_{k}\otimes u_{k})|\mathcal{F}_{k}\Big{]},$ $\displaystyle\leq c_{2}\nabla
f(x_{k})+\mathbb{E}\Big{[}\frac{\delta_{k}^{2}(d+1)u_{k}}{12}(\nabla^{3}f(\bar{x_{k}}^{+})+\nabla^{3}f(\bar{x_{k}}^{-}))(u_{k}\otimes
u_{k}\otimes u_{k})|\mathcal{F}_{k}\Big{]}.$

Now the $j$th coordinate of the second term in RHS of the above inequality is
bounded as follows:

$\displaystyle\mathbb{E}\Big{[}\frac{\delta_{k}^{2}u_{k}^{j}(d+1)}{12}(\nabla^{3}f(\bar{x_{k}}^{+})+\nabla^{3}f(\bar{x_{k}}^{-}))(u_{k}\otimes
u_{k}\otimes u_{k})|\mathcal{F}_{k}\Big{]}$
$\displaystyle\leq\frac{B_{1}\delta_{k}^{2}(d+1)}{6}\sum_{l_{1}=1}^{d}\sum_{l_{2}=1}^{d}\sum_{l_{3}=1}^{d}\mathbb{E}(u_{k}^{j}u_{k}^{l_{1}}u_{k}^{l_{2}}u_{k}^{l_{3}})\leq\frac{B_{1}d^{4}\delta_{k}^{2}}{3}.$

The first inequality follows from (A2) and in the last one we use the fact
$|u_{k}^{l}|\leq 1$. ∎

#### Proof of Theorem 4

###### Proof.

Notice by Lemma 7,

$\mathbb{E}_{\alpha_{[k]}}[\Gamma^{\prime}_{k}]=\mathbb{E}_{\alpha_{[k]}}[\Gamma^{\prime}_{k}|\alpha_{k-1}]=\mathbb{E}_{\alpha_{[k]}}[\Gamma^{\prime}_{k}|x_{k}]=\mathbb{E}_{\alpha[k]}[\tilde{G}-c_{2}\nabla
f|x_{k}]\stackrel{{\scriptstyle(e)}}{{\leq}}\tau\mathbf{1}_{d\times 1},$

with
$\Gamma^{\prime}_{k}\equiv\tilde{G}(x_{k},\xi_{k},\xi_{k}^{+},u_{k},\delta)-c_{2}\nabla
f(x)\equiv\tilde{G}(x_{k},\alpha_{k})-c_{2}\nabla f(x_{k})$ and
$\tau=c_{2}^{\prime\prime}\delta^{2}$ and
$\alpha_{k}\equiv(\xi_{k},\xi_{k}^{+},u_{k},\delta),k\geq 1$. In the above
vector inequality $(e)$ implies that $x_{i}\geq y_{i}$ for
$X,Y\in\mathbb{R}^{d}$. Also, by Section 4.1 we have

$\mathbb{E}_{\alpha_{[k]}}[\lVert\tilde{G}\rVert^{2}]\leq\lVert\mathbb{E}_{\alpha_{[k]}}[\tilde{G}]\rVert^{2}+\frac{C^{\prime\prime}}{\delta^{2}}.$

The rest proof is similar as Theorem 3. ∎

####  4.3.2 Proof of Theorem 6

######  Lemma 8.

Let $\nabla f_{\delta}(x)$ is the balanced estimator defined in 40 then under
(A2) we have

$\lVert\nabla f_{\delta}(x)\rVert^{2}\leq 2\lVert\nabla
f(x)\rVert^{2}c_{11}^{2}+2d^{2}\delta^{4}c_{2}^{\prime\prime 2}$ (43)

###### Proof.

From the definition of balanced estimator we have

$\nabla
f_{\delta}(x_{k})={\displaystyle\mathbb{E}_{u}\bigg{[}\bigg{(}\frac{f(x_{k}+\delta_{k}u_{k})-f(x_{k}-\delta_{k}u_{k})}{2\delta_{k}}\bigg{)}\frac{(d+1)u_{k}}{1+\lVert
u_{k}\rVert^{2}}}\bigg{]}.$

By Taylor series expansion we obtain

$\begin{split}\nabla
f_{\delta}(x_{k})&=\mathbb{E}_{u}\bigg{[}\bigg{(}\frac{2\delta u_{k}^{T}\nabla
f(x_{k})+\frac{\delta^{3}}{6}(\nabla^{3}f(\bar{x_{k}}^{+})+\nabla^{3}f(\bar{x_{k}}^{-}))(u_{k}\otimes
u_{k}\otimes u_{k})}{2\delta}\bigg{)}\frac{(d+1)u_{k}}{1+\lVert
u_{k}\rVert^{2}}\bigg{]}\\\
&=\mathbb{E}_{u}\Big{[}\frac{(d+1)u_{k}u_{k}^{T}\nabla f(x_{k})}{1+\lVert
u_{k}\rVert^{2}}\Big{]}+\mathbb{E}_{u}\Big{[}\frac{\delta^{2}(d+1)u_{k}}{12(1+\lVert
u_{k}\rVert^{2})}(\nabla^{3}f(\bar{x_{k}}^{+})+\nabla^{3}f(\bar{x_{k}}^{-}))(u_{k}\otimes
u_{k}\otimes u_{k})\Big{]},\\\
&\leq\mathbb{E}_{u}\Big{[}\frac{(d+1)u_{k}u_{k}^{T}\nabla f(x_{k})}{1+\lVert
u_{k}\rVert^{2}}\Big{]}+c_{2}^{\prime\prime}\delta^{2}\mathbf{1}_{d}\end{split}$

The second term in the above inequality is obtain via same logic used in Lemma
7. Hence by taking norm of the both side and applying Jensen inequality we get

$\begin{split}\lVert\nabla f_{\delta}(x)\rVert&\leq\lVert\nabla
f(x)\rVert\mathbb{E}\bigg{[}\frac{\lVert u\rVert^{2}(d+1)}{1+\lVert
u\rVert^{2}}\bigg{]}+dc_{2}^{\prime\prime}\delta^{2},\\\ &<\lVert\nabla
f(x)\rVert\mathbb{E}\bigg{[}\lVert
u\rVert^{2}(d+1)\bigg{]}+dc_{2}^{\prime\prime}\delta^{2},\\\ &\leq\lVert\nabla
f(x)\rVert c_{11}+dc_{2}^{\prime\prime}\delta^{2},\\\ \lVert\nabla
f_{\delta}(x)\rVert^{2}&\leq 2\lVert\nabla
f(x)\rVert^{2}c_{11}^{2}+2d^{2}\delta^{4}c_{2}^{\prime\prime 2}.\\\
\end{split}$

∎

#### Proof of Theorem 6

###### Proof.

Define $\alpha_{k}\equiv(\xi_{k},\xi_{k}^{+},u_{k}),k\geq 1$ and
$\alpha_{[N]}:=(\alpha_{1},\alpha_{2},...\alpha_{N})$ as before in Lemma 6
Using Taylor series expansion over $f(x_{k})$ we obtain for any $k=1,2,...,N$,

$\begin{split}f(x_{k+1})&\leq f(x_{k})-\gamma_{k}\langle\nabla
f(x_{k}),\tilde{G}(x_{k},\alpha)\rangle+\frac{L}{2}\gamma_{k}^{2}\left\lVert\tilde{G}(x_{k},\alpha_{k})\right\rVert^{2}\\\
&=f(x_{k})-c_{2}\gamma_{k}\left\lVert\nabla
f(x_{k})\right\rVert^{2}-\gamma_{k}\langle\nabla
f(x_{k}),\Gamma^{\prime}_{k}\rangle+\frac{L}{2}\gamma_{k}^{2}\left\lVert\tilde{G}(x_{k},\alpha_{k})\right\rVert^{2}.\end{split}$

Here
$\Gamma^{\prime}_{k}\equiv\tilde{G}(x_{k},\xi_{k},\xi_{k}^{+},u_{k},\delta)-c_{2}\nabla
f(x)\equiv\tilde{G}(x_{k},\alpha_{k})-c_{2}\nabla f(x_{k})$. Adding upto
$N$-terms both side of these inequalities and applying $f^{*}\leq f(x_{N+1})$,
we obtain

$\displaystyle\sum_{k=1}^{N}c_{2}\gamma_{k}\left\lVert\nabla
f(x_{k})\right\rVert^{2}$ $\displaystyle\leq
f(x_{1})-f^{*}-\sum_{k=1}^{N}\gamma_{k}\langle\nabla
f(x_{k}),\Gamma^{\prime}_{k}\rangle+\frac{L}{2}\sum_{k=1}^{N}\gamma_{k}^{2}\left\lVert\tilde{G}(x_{k},\alpha_{k})\right\rVert^{2}.$

Notice by Lemma 7

$\mathbb{E}_{\alpha_{[k]}}[\Gamma^{\prime}_{k}]=\mathbb{E}_{\alpha_{[k]}}[\Gamma^{\prime}_{k}|\alpha_{k-1}]=\mathbb{E}_{\alpha_{[k]}}[\Gamma^{\prime}_{k}|x_{k}]=\mathbb{E}_{\alpha[k]}[\tilde{G}-c_{2}\nabla
f|x_{k}]\stackrel{{\scriptstyle(e)}}{{\leq}}\tau\mathbf{1}_{d\times 1},$

where $\Gamma^{\prime}_{k}\equiv\tilde{G}(x_{k},\xi_{k},u_{k})-c_{2}\nabla
f(x)\equiv\tilde{G}(x_{k},\alpha_{k})-c_{2}\nabla f(x_{k})$ and
$\tau=c_{2}^{\prime\prime}\delta^{2}$. In the above vector inequality $(e)$
implies that $x_{i}\geq y_{i}$ for $X,Y\in\mathbb{R}^{d}$.

Now by (A9) and Lemma 8

$\begin{split}\mathbb{E}[\left\lVert
G(x_{k},\alpha_{k})\right\rVert^{2}|\alpha_{[k-1]}]&\leq\left[2c^{2}_{11}\mathbb{E}[\lVert\nabla
F(x_{k},\xi_{k})\rVert^{2}|\alpha_{[k-1]}]+2d^{2}\delta^{4}c_{2}^{\prime\prime
2}\right]\\\ &\leq\left[4c^{2}_{11}\mathbb{E}[\left\lVert\nabla
f(x_{k})\right\rVert^{2}|\alpha_{[k-1]}+\sigma^{2}]+2d^{2}\delta^{4}c_{2}^{\prime\prime
2}\right].\end{split}$ (44)

Thus we have

$\displaystyle\sum_{k=1}^{N}c_{2}\gamma_{k}\mathbb{E}_{\alpha_{[N]}}\left\lVert\nabla
f(x_{k})\right\rVert^{2}$ $\displaystyle\leq
D+B\tau\sum_{k=1}^{N}\gamma_{k}+\frac{L}{2}\sum_{k=1}^{N}\gamma_{k}^{2}\left[4c^{2}_{11}\mathbb{E}[\left\lVert\nabla
f(x_{k})\right\rVert^{2}|\alpha_{[k-1]}+\sigma^{2}]+2d^{2}\delta^{4}c_{2}^{\prime\prime
2}\right].$

The above inequality uses the fact $-\lVert
V\rVert_{1}\leq\sum_{k=1}^{d}v_{k}$ for a $d-$dimensional vector $V$ followed
by $\lVert\nabla f(x_{k})\rVert_{1}\leq\lVert\nabla f(x_{k})\rVert\leq B$ in
(A7). Note that $D=f(x_{1})-f^{*}$. By rearranging the terms we have

$\displaystyle\sum_{k=1}^{N}\left[c_{2}\gamma_{k}-2Lc^{2}_{11}\gamma_{k}^{2}\right]\mathbb{E}_{\alpha_{[N]}}\left\lVert\nabla
f(x_{k})\right\rVert^{2}$ $\displaystyle\leq
D+B\tau\sum_{k=1}^{N}\gamma_{k}+L\left(4c_{11}\sigma^{2}+2d^{2}\delta^{4}c_{2}^{\prime\prime
2}\right)\sum_{k=1}^{N}\gamma_{k}^{2}.$

By the same argument as in 35 and under the probability distribution

$\mathbb{P}_{R}(k)=Prob(R=k)=\frac{\left[c_{2}\gamma_{k}-2Lc^{2}_{11}\gamma_{k}^{2}\right]}{\sum_{k=1}^{N}\left[c_{2}\gamma_{k}-2Lc^{2}_{11}\gamma_{k}^{2}\right]},$
(45)

the following is obtained

$\displaystyle\mathbb{E}_{\alpha_{[N]}}\left\lVert\nabla
f(x_{R})\right\rVert^{2}$
$\displaystyle\leq\frac{1}{\sum_{k=1}^{N}\left[c_{2}\gamma_{k}-Lc^{2}_{11}\gamma_{k}^{2}\right]}\bigg{[}D+B\tau\sum_{k=1}^{N}\gamma_{k}+L\left(4c_{11}\sigma^{2}+2d^{2}\delta^{4}c_{2}^{\prime\prime
2}\right)\sum_{k=1}^{N}\gamma_{k}^{2}.\bigg{]}.$

Note that, the condition of $\gamma_{k}$, see 15, is given by

$\gamma_{k}=\min\left\\{\frac{c_{2}}{4c^{2}_{11}L},\frac{1}{N^{1/2}}\right\\},\quad
k=1,2,...N,$

Thus one can have

$\sum_{k=1}^{N}\left[c_{2}\gamma_{k}-2Lc^{2}_{11}\gamma_{k}^{2}\right]=Nc_{2}\gamma_{1}\left[1-\frac{2Lc^{2}_{11}}{c_{2}}\gamma_{1}\right]\geq\frac{Nc_{2}\gamma_{1}}{2}.$

From the above inequality, we can write

$\displaystyle\mathbb{E}\left[\left\lVert\nabla
f(x_{R})\right\rVert^{2}\right]$
$\displaystyle\leq\frac{2D}{Nc_{2}\gamma_{1}}+\frac{2B\tau}{c_{2}}+\frac{L}{c_{2}}\left(4c_{11}\sigma^{2}+2d^{2}\delta^{4}c_{2}^{\prime\prime
2}\right)\gamma_{1},$
$\displaystyle\leq\frac{2D}{Nc_{2}}\max\left\\{\frac{4c^{2}_{11}L}{c_{2}},N^{1/2}\right\\}+\frac{2B\tau}{c_{2}}+\frac{L}{c_{2}N^{1/2}}\left(4c_{11}\sigma^{2}+2d^{2}\delta^{4}c_{2}^{\prime\prime
2}\right),$
$\displaystyle\stackrel{{\scriptstyle(h)}}{{\leq}}\Big{(}\frac{2DL}{Nc_{2}^{2}}+\frac{2D}{N^{1/2}}\Big{)}+\frac{2Bc_{2}^{\prime\prime}}{c_{2}N}+\frac{L}{2c_{2}N^{1/2}}\left(4c_{11}\sigma^{2}+\frac{2d^{2}c_{2}^{\prime\prime}}{N^{2}}\right).$

Note that $(h)$ uses the condition of $\delta=\frac{1}{N^{1/2}}$ by putting
$\tau=c_{2}^{\prime\prime}\delta^{2}$. Thus, we get a rate of convergence of
$\mathcal{O}(N^{-1/2}).$ ∎

##  5 Experiments

In this section, we compare the performance of GSF, SPSA with symmetric
Bernoulli($\pm 1$) valued perturbations and RDSA with uniform ($-5,5$)
perturbations, respectively, with the TCSF (Algorithm 1) and TCSF with
balanced estimator13 (B-TCSF). We consider both non-convex and convex
objective functions with additive noise. We consider the following choices for
the function $F(x,\xi)$ with $d=4$:

Table 1: The different models for the observations $F(\cdot,\cdot)$ considered in our experiments. Name | Functional form | Optimal point $x^{*}$ | Min-Value
---|---|---|---
Rastrigin | $10d+\sum_{i=1}^{d}(x_{i}^{2}-10\cos(2\pi x_{i}))+\xi_{x}$ | $(0,\ldots,0)^{T}$ | 0
Rosenbrock |  $\sum_{i=1}^{d-1}(100(x_{i+1}-x_{i}^{2})^{2}+(1-x_{i})^{2})+\xi_{x}$, | $(1,\ldots,1)^{T}$ | 0
Quadratic | $\frac{1}{2}x^{T}Ax-b^{T}x+\xi_{x}$ | $(-135.1,\ldots,-5.6)^{T}$ | $-17.528$

The setting for the case of quadratic function considered in Table 1 is

$F_{3}(x,\xi_{x})=\frac{1}{2}x^{T}Ax-b^{T}x+\xi_{x},$

with, $x^{*}=[-135.1150,-4.5224,130.1168,-5.6879]^{T}$ where

$A=\begin{bmatrix}2.3346&1.1384&2.5606&14507\\\ 1.1384&0.7860&1.2743&0.9531\\\
2.5606&1.2743&2.8147&1.6487\\\ 1.4507&0.9531&1.6487&1.8123\end{bmatrix}$

and

$b=[0.4218,0.9157,0.7922,0.9595]^{T}.$

We consider three settings for the noise $\eta$. In the first setting,
referred to as Type-1, we let $\xi_{x}=[x^{T},1]\eta$, where $\eta$ is a
multivariate Gaussian with zero mean, and covariance matrix
$\sigma^{2}\mathbb{I}_{d+1}$ with $\sigma=5$. In the second setting, referred
to as Type-2, we have $\xi_{x}$ as a Gaussian random variable with mean $0$
and variance $\ln\lVert x\rVert_{2}$. Finally, in the last setting, referred
to as Type-3, we have $\xi_{x}$ as a Gaussian random variable with mean $0$
and variance $\frac{1}{1+\ln\lVert x\rVert_{2}}$. Note in particular that
Rosenbrock is a badly-scaled function, while Rastrigin is a multi-modal
function.

For the truncated Cauchy perturbations in TCSF and B-TCSF, we generated
samples from the multivariate t-distribution with one degree of freedom and
then projected the same to the unit sphere. In our experiments, we set
$\epsilon=0.0001$, i.e., we stop the algorithm when $\lVert
G\rVert\leq\epsilon$. For our initial experiments, we used the stepsize and
smoothing parameter as follows: $\gamma_{k}=\frac{1}{k^{0.6}}$ and
$\delta_{k}=\frac{1}{k^{0.09}}$. However, we used constant step-sizes $0.0001$
and smoothing parameters $0.001$, respectively, which work well on both
algorithms. We run each algorithm 100 times with $N=1000$, $3000$, $10000$,
respectively, for the Rastrigin, quadratic and Rosenbrock functions, and the
averages of the optimal functional values are reported in Section 5 while the
standard error estimates for the various algorithms from the different
simulation runs are given in Section 5. we considered the set $[0,10]^{4}$ for
Rastrigin and Rosenbrock functions, and $[0,150]^{4}$ for a quadratic
objective function. We report the average number of iterations needed to reach
an $\epsilon$-stationary point in Section 5

Table 2: Average functional values for five SG algorithms with diminishing
step size and smoothing parameter

|  |  GSF |  TCSF |  B-TCSF |  SPSA |  RDSA
---|---|---|---|---|---|---
| Rastrigin | $0.00019$ | $1.03e-05$ | $0.0$ | $0.0092$ | $0.0094$
| Rosenbrock | $0.002$ | $0.00066$ | $0.00062$ | $0.0010$ | $0.0017$
Type 1 | Quadratic | $-17.471$ | $-17.513$ | $-17.5286$ | $-17.4719$ | $-17.4731$
| Rastrigin | $0.0097$ | $0.0009$ | $1.17e-05$ | $0.0096$ | $0.01$
| Rosenbrock | $0.003$ | $0.0018$ | $0.00091$ | $8e-05$ | $0.0017$
Type 2 | Quadratic | $-17.5286$ | $-17.528$ | $-17.5248$ | $-17.52819$ | $-17.474$
| Rastrigin | $0.0002$ | $1.11e-5$ | $0.0$ | $0.0$ | $0.009$
| Rosenbrock | $0.017$ | $0.00072$ | $0.00026$ | $0.0021$ | $0.0017$
Type 3 | Quadratic | $-17.5253$ | $-17.5248$ | $-17.5279$ | $-17.492$ | $-17.488$
|  |  |  |  |  |

Table 3: Standard error for Section 5

|  |  RSGF |  TCSF |  B-TCSF |  SPSA |  RDSA
---|---|---|---|---|---|---
| Rastrigin | $3.44\times 10^{-7}$ | $2.34\times 10^{-8}$ | $0$ | $9.41\times 10^{-6}$ | $4.84\times 10^{-6}$
| Rosenbrock | $4.96\times 10^{-5}$ | $3.75\times 10^{-7}$ | $1.24\times 10^{-7}$ | $6.34\times 10^{-6}$ | $1.51\times 10^{-6}$
Type 1 | Quadratic | $5.61\times 10^{-6}$ | $1.58\times 10^{-7}$ | $1.28\times 10^{-7}$ | $64.25\times 10^{-5}$ | $5.28\times 10^{-6}$
| Rastrigin | $1.26\times 10^{-5}$ | $1.62\times 10^{-7}$ | $2.96\times 10^{-7}$ | $4.29\times 10^{-5}$ | $5.61\times 10^{-6}$
| Rosenbrock | $1.59\times 10^{-6}$ | $4.57\times 10^{-7}$ | $4.29\times 10^{-7}$ | $4.27\times 10^{-6}$ | $1.59\times 10^{-6}$
Type 2 | Quadratic | $4.36\times 10^{-7}$ | $1.54\times 10^{-7}$ | $2.03\times 10^{-6}$ | $6.29\times 10^{-6}$ | $5.69\times 10^{-6}$
| Rastrigin | $4.29\times 10^{-6}$ | $1.80\times 10^{-7}$ | $0$ | $0$ | $2.49\times 10^{-5}$
| Rosenbrock | $2.41\times 10^{-6}$ | $1.94\times 10^{-5}$ | $1.62\times 10^{-7}$ | $2.36\times 10^{-5}$ | $6.16\times 10^{-5}$
Type 3 | Quadratic | $3.26\times 10^{-6}$ | $2.63\times 10^{-6}$ | $1.30\times 10^{-7}$ | $6.11\times 10^{-5}$ | $8.83\times 10^{-5}$
|  |  |  |  |  |

There is a significant difference in optimal functional values obtained from
TCSF, B-TCSF as compared to the other algorithms. One can notice from the
Section 5 that $|f(x^{*})-f(\bar{x}_{sol})|$ $>$
$|f(x^{*})-f(\hat{x}_{{sol}})|$, where $\bar{x}_{sol}$ indicates the final
output from GSF, SPSA, RDSA and $\hat{x}_{{sol}}$ denotes the final output
from TCSF, B-TCSF. Only in the case of Rosenbrock and Rastrigin with Type 2
and Type 3 error, SPSA works slightly better than TCSF though not so when
compared with B-TCSF. However in case of the quadratic function we have
noticed that TCSF beats other algorithms under Type-1 and Type-3 errors as
well. Between TCSF and B-TCSF, B-TCSF is seen to perform slightly better on
the whole. This is to be expected since B-TCSF has a lower bias because of the
direct cancellation of even-order bias terms starting from the second order.

Table 4: Average number of iterations to converge to the optimal
point(constant step size)

|  |  RSGF |  TCSF |  B-TCSF |  SPSA |  RDSA
---|---|---|---|---|---|---
| Rastrigin | $837.6$ | $331.5$ | $149.6$ | $1268.39$ | $1281.39$
| Rosenbrock | $5604.08$ | $3995.2$ | $3138.3$ | $6968.41$ | $8572.53$
Type 1 | Quadratic | $8727.82$ | $3234.86$ | $2994.15$ | $5952.7$ | $6000.3$
| Rastrigin | $785.6$ | $638.9$ | $151.4$ | $1269.32$ | $1105.33$
| Rosenbrock | $5924.5$ | $4834.2$ | $2231.2$ | $1782.1$ | $9419.4$
Type 2 | Quadratic | $3305.44$ | $3050.96$ | $2863.11$ | $3154.49$ | $5743.9$
| Rastrigin | $896.8$ | $398.6$ | $128.7$ | $137.6$ | $1352.4$
| Rosenbrock | $5450.51$ | $3376.6$ | $2735.2$ | $7946.4$ | $8847.3$
Type 3 | Quadratic | $3518.9$ | $3319.68$ | $2726.65$ | $4154.49$ | $6152.2$
|  |  |  |  |  |

Table 5: Standard error for Section 5

|  |  RSGF |  TCSF |  B-TCSF |  SPSA |  RDSA
---|---|---|---|---|---|---
| Rastrigin | $4.79$ | $0.325$ | $0.1255$ | $0.395$ | $0.445$
| Rosenbrock | $90.4$ | $38.35$ | $20.3$ | $108.2$ | $115.3$
Type 1 | Quadratic | $31.85$ | $24.49$ | $17.85$ | $1.07$ | $1.0407$
| Rastrigin | $5.66$ | $5.412$ | $0.14$ | $0.445$ | $0.415$
| Rosenbrock | $108.08$ | $45.36$ | $21.71$ | $13.99$ | $99.35$
Type 2 | Quadratic | $20.78$ | $20.85`$ | $15.18$ | $0.895$ | $1.13$
| Rastrigin | $6.114$ | $0.319$ | $0.132$ | $0.128$ | $0.431$
| Rosenbrock | $102.21$ | $22.93$ | $19.02$ | $103.36$ | $121.3$
Type 3 | Quadratic | $22.17$ | $24.18$ | $19.15$ | $1.04$ | $0.99$
|  |  |  |  |  |

In Section 5 above we have described the number of iterations needed to reach
an $\epsilon-$stationary point while Section 5 describes the SE for
corresponding iteration . There is a significant difference in the number of
iterations required for converging to the optimal point in each of the cases.
For example, GSF, RDSA and SPSA each take more than 2000 iterations with
respect to TCSF, B-TCSF. However SPSA performs well as compared to RDSA, RSGF
and in some experiments it beats TCSF (e.g., Rosenbrock with type-2 error).

Table 6: Average functional values for five SG algorithm with constant step
size and smoothing parameter

|  |  RSGF |  TCSF |  B-TCSF |  SPSA |  RDSA
---|---|---|---|---|---|---
| Rastrigin | $1.722$ | $0.091$ | $0.008$ | $0.987$ | $2.024$
| Rosenbrock | $1.954$ | $0.013$ | $0.0006$ | $0.587$ | $1.743$
Type 1 | Quadratic | $-15.894$ | $-17.685$ | $-17.587$ | $-17.096$ | $-14.843$
| Rastrigin | $1.985$ | $0.065$ | $0.008$ | $0.848$ | $2.542$
| Rosenbrock | $1.926$ | $0.0018$ | $0.0089$ | $0.683$ | $2.193$
Type 2 | Quadratic | $-14.448$ | $-17.939$ | $-17.373$ | $-16.597$ | $-13.2493$
| Rastrigin | $1.875$ | $0.01$ | $0.00623$ | $0.879$ | $1.893$
| Rosenbrock | $2.097$ | $0.091$ | $0.0034$ | $0.698$ | $2.314$
Type 3 | Quadratic | $-14.869$ | $-17.962$ | $-17.459$ | $-16.183$ | $-13.743$
|  |  |  |  |  |

We can observe from Section 5 that for constant step size and smoothing
parameter $f(\bar{x}_{TCSF})$,$f(\bar{x}_{B-TCSF})$ provides almost optimal
functional value as compared to other algorithms. However, B-TCSF is more
efficient than TCSF in this context. Thus we can conclude that there is an
improvement in empirical performance when TCSF and B-TCSF are used over other
SG algorithms.

##  6 Conclusions and future work

We proposed and analyzed a gradient estimation scheme, based on truncated
Cauchy random perturbations, for solving a non-convex smooth optimization
problem. We showed that our algorithm avoids traps and converges
asymptotically to a local minimum. Our algorithm performs better than two
popular gradient estimation schemes in the literature, namely SPSA and GSF, in
terms of the asymptotic convergence rate. We also provided non-asymptotic rate
for our algorithm that is the same as the asymptotic rate and better when
common random noise is used in the simulations. Our algorithm also performs
better than GSF, SPSA, RDSA empirically. Exploring the performance of the
Newton method using Hessian estimation under the truncated Cauchy
perturbations would be an interesting direction for future work.

## References

  * Baes (2009) M. Baes.  Estimate sequence methods: Extensions and approximations.  _IFOR Internal report, ETH Zurich, Switzerland, Tech_ , 2009.
  * Balasubramanian and Ghadimi (2018) K. Balasubramanian and S. Ghadimi.  Zeroth-order nonconvex stochastic optimization: Handling constraints, high dimensionality, and saddle points.  _Springer_ , 2018.
  * Bertsekas (2019) D. Bertsekas.  _Reinforcement learning and optimal control_.  Athena Scientific, 2019.
  * Bhatnagar (2007) S. Bhatnagar.  Adaptive Newton-based smoothed functional algorithms for simulation optimization.  _ACM Transactions on Modeling and Computer Simulation_ , 18(1):2:1–2:35, 2007.
  * Bhatnagar and Borkar (2003) S. Bhatnagar and V. S. Borkar.  Multiscale chaotic spsa and smoothed functional algorithms for simulation optimization.  _Simulation 79_ , page 568–580, 2003.
  * Bhatnagar et al. (2013) S. Bhatnagar, H. L. Prasad, and L. A. Prashanth.  _Stochastic Recursive Algorithms for Optimization: Simultaneous Perturbation Methods (Lecture Notes in Control and Information Sciences)_.  Springer, 2013.
  * Bhavsar and Prashanth (2022) N. Bhavsar and L. A. Prashanth.  Non-asymptotic bounds for stochastic optimization with biased noisy gradient oracles.  _IEEE Transactions on Automatic Control_ , 2022.
  * Chin (1997) D. C. Chin.  Comparative study of stochastic algorithms for system optimization based on gradient approximations.  _IEEE Transactions on Systems, Man, and Cybernetics, Part B: Cybernetics_ , 27(2):244–249, 1997.
  * Devolder et al. (2014) O. Devolder, F. Glineur, and Y. Nesterov.  First- order methods of smooth convex optimization with inexact oracle.  _Mathematical Programming_ , page 37–75, 2014.
  * Gawthrop and Sbarbaro (1990) P. J. Gawthrop and D. Sbarbaro.  Stochastic approximation and multilayer perceptrons: The gain backpropagation algorithm.  _Complex Systems_ , 4(1):51–74, 1990.
  * Ge et al. (2015) R. Ge, F. Huang, C. Jin, and Y. Yuan.  Escaping from saddle points – online stochastic gradient for tensor decomposition.  _JMLR: Workshop and Conference Proceedings_ , 2015.
  * Gerencser (1999) L. Gerencser.  Convergence rate of moments in stochastic approximation with simultaneous perturbation gradient approximation and resetting.  _IEEE Trans. Autom. Control_ , page 894–905, 1999.
  * Ghadimi and Lan (2013) S. Ghadimi and G. Lan.  Stochastic first- and zeroth-order methods for nonconvex stochastic programming.  _SIAM Journal on Optimization_ , pages 15–20, 2013.
  * Hu et al. (2016) X. Hu, L. A. Prashanth, A. Gyorgy, and C. Szepesvari.  (bandit) convex optimization with biased noisy gradient oracles.  _Artificial Intel- ligence and Statistics_ , page 819–828, 2016\.
  * Jin et al. (2017) C. Jin, R. Ge, P. Netrapalli, S. M. Kakade, and M. I. Jordan.  How to escape saddle points efficiently.  _ICML_ , 2017.
  * Katkovnik and Kulchitsky (1972) V. Y. Katkovnik and Y. Kulchitsky.  Convergence of a class of random search algorithms.  _Automation Remote Control_ , 8:1321–1326, 1972.
  * Kiefer and Wolfowitz (1952) J. Kiefer and J. Wolfowitz.  Stochastic estimation of the maximum of a regression function.  _The Annals of Mathematical Statistics_ , page 462–466, 1952.
  * Kreimer and Rubinstein (1972) J. Kreimer and R. Y. Rubinstein.  Smoothed functionals and constrained stochastic approximation.  _JSTOR,_ , 1972.
  * Kushner and Clark (1978a) H. Kushner and D. Clark.  Stochastic approximation methods for constrained and unconstrained systems.  _Springer Verlag_ , 1978a.
  * Kushner and Clark (1978b) H. J. Kushner and D. S. Clark.  _Stochastic Approximation Methods for Constrained and Unconstrained Systems_.  Springer Verlag, New York, 1978b.
  * Nesterov and Spokoiny (2017) Y. Nesterov and V. Spokoiny.  Random gradient-free minimization of convex functions.  _Foundations of Computational Mathematics volume_ , page 527–566, 2017.
  * Pemantle (1990) R. Pemantle.  Non-convergence to unstable points in urn models and stochastic approximations.  _The annals of Probability_ , page 3, 1990.
  * Prashanth et al. (2017) L. Prashanth, S. Bhatnagar, M. C. Fu, and S. I. Marcus.  Adaptive system optimization using random directions stochastic approximation.  _IEEE Transactions on Automatic Control_ , 62(5):2223–2238, 2017.
  * Robbins and Monro (1951) H. Robbins and S. Monro.  A stochastic approximation method.  _Annals of Mathematical Statistics_ , pages 400–407, 1951.
  * Rubinstein (1981) R. Y. Rubinstein.  _Simulation and the Monte Carlo Method_.  Wiley, New York, 1981.
  * Spall (1992) J. C. Spall.  Multivariate stochastic approximation using a simultaneous perturbation gradient approximation.  _IEEE TRANSACTIONS ON AUTOMATIC CONTROL,_ , 1992.
  * Staneski (1990) P. G. Staneski.  The truncated cauchy distribution: Estimation of parameters and application to stock returns.  _Dissertation. Old Dominion University_ , 1990.
  * Styblinski and Tang (1990) M. Styblinski and T. Tang.  Experiments in non-convex optimization: Stochastic approximation with function smoothing and simulated annealing.  _Neural Networks_ , 1990.
  * Yao et al. (2009) H. Yao, S. Bhatnagar, and C. Szepesvári.  Lms-2: Towards an algorithm that is as cheap as LMS and almost as efficient as RLSf.  In _Proceedings of the 48h IEEE Conference on Decision and Control (CDC) held jointly with 2009 28th Chinese Control Conference_ , pages 1181–1188. IEEE, 2009.
  * Zhu (2021) J. Zhu.  Hessian estimation via stein’s identity in black-box problems.  _Proceedings of Machine Learning Research 2021_ , pages 1–17, 2021\.