# Theoretical Study and Comparison of SPSA and RDSA Algorithms with Different
Perturbations

Yiwen Chen

###### Abstract

Stochastic approximation (SA) algorithms are widely used in system
optimization problems when only noisy measurements of the system are
available. This paper studies two types of SA algorithms in a multivariate
Kiefer-Wolfowitz setting: random-direction SA (RDSA) and simultaneous-
perturbation SA (SPSA), and then describes the bias term, convergence, and
asymptotic normality of RDSA algorithms. The gradient estimations in RDSA and
SPSA have different forms and, consequently, use different types of random
perturbations. This paper looks at various valid distributions for
perturbations in RDSA and SPSA and then compares the two algorithms using
mean-square errors computed from asymptotic distribution. From both a
theoretical and numerical point of view, we find that SPSA generally
outperforms RDSA.

††1\. A compressed version of this paper appears in the 2021 55rd Annual
Conference on Information Sciences and Systems (CISS), Mar 2021.
2\. Contact Information<EMAIL_ADDRESS>/<EMAIL_ADDRESS>
##  1 Introduction

Stochastic approximation is a well-known recursive procedure for finding roots
of equations in the presence of noisy measurements. Consider a smooth
$p$-dimensional loss function $L:\,\mathbb{R}^{p}\rightarrow\mathbb{R}$, with
gradient $\bm{g}:\mathbb{R}^{p}\rightarrow\mathbb{R}^{p}$. Assume that $L$ has
a unique minimum $\bm{\uptheta}^{*}\in\mathbb{R}^{p}$; i.e.
$L(\bm{\uptheta}^{*})\leq L(\bm{\uptheta})$ for all
$\bm{\uptheta}^{*}\in\mathbb{R}^{p}$, and $\bm{g}(\bm{\uptheta}^{*})=\bm{0}$.
In many cases where only noisy measurements of the gradient are available, the
Robbins-Monro stochastic approximation (SA) algorithm is widely used with the
form of:

$\displaystyle\bm{\hat{\uptheta}}_{k+1}=\bm{\hat{\uptheta}}_{k}-a_{k}\bm{Y}_{k}(\bm{\hat{\uptheta}}_{k}).\;k=0,1,2,...$
(1)

where
$\bm{Y}_{k}(\bm{\hat{\uptheta}}_{k})=\partial{Q}/\partial{\bm{\uptheta}}$
represents a direct noisy measurement of the true gradient
$\bm{g}_{k}(\bm{\hat{\uptheta}}_{k})$ for $Q$ in the representation
$L(\bm{\uptheta})=E[Q(\bm{\uptheta,V})]$ with $\bm{V}$ corresponding to the
randomness in the problem. $a_{k}>0$ is the step size, determining the
convergence of the algorithm.

However, in most practical cases, such direct measurement is difficult to
obtain and only noisy measurements of the loss function are available. Then it
comes to the use of gradient-free algorithms. There are several common
algorithms such as the finite-difference (FD) SA, random-direction SA (RDSA)
and the simultaneous-perturbation SA (SPSA). All of them use the measurement
of the loss function to estimate the gradient. The FDSA algorithm has also
been discussed comprehensively in [15]. It has an $i$th component in the form
of:

$\displaystyle\hat{\bm{g}}_{ki}(\hat{\bm{\uptheta}}_{k})=\dfrac{y(\hat{\bm{\uptheta}}_{k}+c_{k}\bm{e}_{i})-y(\hat{\bm{\uptheta}}_{k}-c_{k}\bm{e}_{i})}{2c_{k}},$
(2)

where $\bm{e}_{i}$ denotes the unit vector along the $i$th axis and $y$ is the
noisy measurement of the loss value and $c_{k}>0$ defines the difference
magnitude. The pair ${\\{a_{k},c_{k}}\\}$ represents the gain sequences.
Nevertheless, the FDSA algorithm shows inefficiency in solving problems with
high dimension due to its use of $2p$ measurements of $L$ per iteration. The
other two RDSA and SPSA approaches, alleviate this problem by requiring only
two system simulations regardless of the parameter dimension.

Several papers have discussed the three SA algorithms. Spall has discussed
details of FDSA and SPSA algorithms in [13] and [15]; Kushner discussed the
RDSA algorithm in [9]; Blakney and Zhu compared FDSA and SPSA algorithms in
[2] and Chin compared the three algorithms using mean-square errors in [3].
Because the FDSA and SPSA algorithms have been well analyzed, this paper will
mainly focus on the RDSA algorithm. Some papers investigated the theoretical
foundation for RDSA algorithms such as [1], [18] and [3]; however, some of
them gave the illustration only briefly. This paper is going to arrange and
organize those different versions of theoretical proofs to show the
convergence and asymptotic normality. After that, we will conduct the
comparison between RDSA and SPSA using Mean-Square Error and present certain
conditions under which SPSA outperforms RDSA with specific perturbations.

The remainder of this paper is organized as follows. Section 2 introduces the
general formulation of the SPSA and RDSA algorithms. Section 3 covers the bias
term, convergence condition, and asymptotic normality of RDSA algorithm and
leads to the comparison of relative accuracy of SPSA and RDSA. Section 4
presents several numerical studies to assist the theoretical conclusion in
Section 3. Especially, this paper uses the skewed-quartic function as an
example to show how the theory works on the performance of distinct
distributed perturbation in SPSA and RDSA.

##  2 Background

This section briefly discusses the SPSA algorithm and the RDSA algorithm with
their estimates for the gradient, $\hat{\bm{g}}(\hat{\bm{\uptheta}}_{k})$.

SPSA (‘Simultaneous Perturbation Stochastic Algorithm’) has been well
discussed in [15, Chap. 7] and [14]. The estimate for
$\bm{g}(\hat{\bm{\uptheta}}_{k})$ is in the form of:

$\displaystyle\hat{\bm{g}}_{k}(\hat{\bm{\uptheta}}_{k})=\dfrac{y(\hat{\bm{\uptheta}}_{k}+c_{k}\bm{\Delta}_{k})-y(\hat{\bm{\uptheta}}_{k}-c_{k}\bm{\Delta}_{k})}{2c_{k}\bm{\Delta}_{k}},$

where $y(\hat{\bm{\uptheta}}_{k}\pm
c_{k}\bm{\Delta}_{k})=L(\hat{\bm{\uptheta}}_{k}\pm
c_{k}\bm{\Delta}_{k})+\upepsilon^{(\pm)}_{k}$ with $\upepsilon^{(\pm)}_{k}$
representing measurement noise term that satisfies
$E[\upepsilon^{(+)}_{k}-\upepsilon^{(-)}_{k}]=0$.
$\bm{\Delta}_{k}\in\mathbb{R}^{p}$ is a vector of $p$ mutually independent
mean-zero random variables representing the perturbation, which satisfies
certain conditions, and $\bm{\Delta}^{-1}_{k}$ represents the vector of
inverse components. There’s no specific type of distribution for
$\bm{\Delta}_{k}$ as long as the perturbation satisfies general conditions
mentioned in [15, Chap. 7] (Principally, mean 0, symmetric, finite moments and
certain inverse moments). Usually, $\bm{\Delta}_{k}$ can be chosen as
Bernoulli, U-shape, and other distributions. The range of choices is
restricted due to the main constraint that the inverse moment
$E(|1/\Delta_{km}|)$ should be finite. Many common distributions, such as
Gaussian distribution and Uniform distribution, cannot be applied to the SPSA
algorithm since they contain too much probability mass near 0.

For RDSA (‘Random Direction Stochastic Algorithm’), the estimate for
$\bm{g}(\hat{\bm{\uptheta}}_{k})$ is in the form of:

$\displaystyle\hat{\bm{g}}_{k}(\hat{\bm{\uptheta}}_{k})=\dfrac{y(\hat{\bm{\uptheta}}_{k}+c_{k}\bm{\pi}_{k})-y(\hat{\bm{\uptheta}}_{k}-c_{k}\bm{\pi}_{k})}{2c_{k}}\bm{\pi}_{k},$

where $\bm{\pi}_{k}\in\mathbb{R}^{p}$ is a vector of Monte-Carlo generated
random variables satisfying certain regularity conditions and normalized so
that $||\bm{\pi}_{k}||^{2}=p$ . Unlike SPSA, the requirement for finite
inverse moment is released for $\bm{\pi}_{k}$ in RDSA, so the range of choices
of distribution for $\bm{\pi}_{k}$ is different from that for
$\bm{\Delta}_{k}$ in SPSA. Several valid choices have been discussed in [18]
such as axis distribution, Gaussian distribution, and Uniformly Spherical
distribution where the perturbation distributes uniformly on a $p$-dimensional
sphere.

The main difference between the two algorithms is that in RDSA, both
adjustments to $\hat{\bm{\uptheta}}_{k}$ and
$\hat{\bm{g}}_{k}(\hat{\bm{\uptheta}}_{k})$ are in the same direction
$\bm{\pi}_{k}$ while the SPSA employs two different directions,
$\bm{\Delta}_{k}$ and $\bm{\upzeta}_{k}$ with component
$\upzeta_{ki}=1/{\Delta}_{ki}$.

There are two popular choices for $\bm{\pi}_{k}$: independent Gaussian
distribution $N(0,1)$ and Uniformly Spherical distribution with radius
$\sqrt{p}$. Notice that in Uniformly Spherical distribution, since there is
the normalization, $||\bm{\pi}_{k}||^{2}=p$, variables are not independent,
but uncorrelated.

##  3 Methodology

This section presents several results that form the theoretical basis for the
RDSA algorithm, similar to those results of SPSA, which have been completely
proved in [12]. The following sections consider the bias term in
$\bm{\hat{g}}_{k}(\bm{\hat{\uptheta}}_{k})$, the strong convergence and
asymptotic normality of $\hat{\bm{\uptheta}}_{k}$ in the RDSA algorithm. Then
we compare the performance of RDSA with SPSA using mean-square error computed
with the asymptotic distribution.

###  3.1 The bias in $\hat{\bm{g}}_{k}$

The bias term in $\hat{\bm{g}}_{k}$ is in the form of:

$\displaystyle\bm{b}_{k}(\hat{\bm{\uptheta}}_{k})=E[\hat{\bm{g}}_{k}(\hat{\bm{\uptheta}}_{k})|\hat{\theta}_{k}]-\bm{g}_{k}(\bm{\hat{\uptheta}}_{k}),$
(3)

where $\upepsilon^{(\pm)}_{k}=y(\bm{\uptheta}_{k}\pm
c_{k}\bm{\pi}_{k})-L(\bm{\uptheta}_{k}\pm c_{k}\bm{\pi}_{k})$ is the noise on
the estimate of the loss function $L$, with mean 0 and variance
$\upsigma^{2}$. Here we henceforth assume that the noise is independent with
$\hat{\bm{\uptheta}}_{k}$ and $\bm{\pi}_{k}$.

Calculate $E[\hat{\bm{g}}_{k}(\hat{\bm{\uptheta}}_{k})]$ and
$\bm{g}(\hat{\bm{\uptheta}}_{k})$ respectively:

$\displaystyle\hat{\bm{g}}_{k}(\hat{\bm{\uptheta}}_{k})=$
$\displaystyle\dfrac{y(\hat{\bm{\uptheta}}_{k}+c_{k}\bm{\pi}_{k})-y(\hat{\bm{\uptheta}}_{k}-c_{k}\bm{\pi}_{k})}{2c_{k}}\;\bm{\pi}_{k}$
$\displaystyle=\dfrac{L(\hat{\bm{\uptheta}}_{k}+c_{k}\bm{\pi}_{k})+\upepsilon^{(+)}_{k}-L(\hat{\bm{\uptheta}}_{k}-c_{k}\bm{\pi}_{k})-\upepsilon^{(-)}_{k}}{2c_{k}}\;\bm{\pi}_{k}$
$\displaystyle=\dfrac{L(\hat{\bm{\uptheta}}_{k}+c_{k}\bm{\pi}_{k})-L(\hat{\bm{\uptheta}}_{k}-c_{k}\bm{\pi}_{k})}{2c_{k}}\;\bm{\pi}_{k}+\dfrac{\upepsilon^{(+)}_{k}-\upepsilon^{(-)}_{k}}{2c_{k}}\bm{\pi}_{k}$
$\displaystyle=\dfrac{2c_{k}\bm{g}(\hat{\bm{\uptheta}}_{k})^{T}\bm{\pi}_{k}+\dfrac{1}{6}c^{3}_{k}L^{{}^{\prime\prime\prime}}(\bar{\bm{\uptheta}}^{(+)}_{k})[\bm{\pi}_{k}\otimes\bm{\pi}_{k}\otimes\bm{\pi}_{k}]}{2c_{k}}\;\bm{\pi}_{k}$
$\displaystyle+\dfrac{\dfrac{1}{6}c^{3}_{k}L^{{}^{\prime\prime\prime}}(\bar{\bm{\uptheta}}^{(-)}_{k})[\bm{\pi}_{k}\otimes\bm{\pi}_{k}\otimes\bm{\pi}_{k}]}{2c_{k}}\;\bm{\pi}_{k}$
$\displaystyle+O(c^{5}_{k})+\dfrac{\upepsilon^{(+)}_{k}-\upepsilon^{(-)}_{k}}{2c_{k}}\bm{\pi}_{k}$

where $\bar{\bm{\uptheta}}^{(\pm)}_{k}$ denotes points on the line segments
between $\hat{\bm{\uptheta}}_{k}$ and $\hat{\bm{\uptheta}}_{k}\pm
c_{k}\bm{\pi}_{k}$.

Then look at the $m$th term of $\hat{\bm{g}}_{k}(\hat{\bm{\uptheta}}_{k})$:

$\displaystyle\hat{g}_{km}(\hat{\bm{\uptheta}}_{k})=$
$\displaystyle\left(\sum_{i=1}^{p}g_{i}(\hat{\bm{\uptheta}}_{k}){\pi}_{ki}+\dfrac{1}{12}c^{2}_{k}[L^{{}^{\prime\prime\prime}}(\bar{\bm{\uptheta}}^{(+)}_{k})+L^{{}^{\prime\prime\prime}}(\bar{\bm{\uptheta}}^{(-)}_{k})][\bm{\pi}_{k}\otimes\bm{\pi}_{k}\otimes\bm{\pi}_{k}]\right)\,\pi_{km}$
$\displaystyle+\dfrac{\upepsilon^{(+)}_{km}-\upepsilon^{(-)}_{km}}{2c_{k}}\pi_{km}$
$\displaystyle=g_{m}(\hat{\uptheta}_{k})\pi^{2}_{km}+\sum_{i=1,i\neq
m}^{p}g_{i}(\hat{\uptheta}_{k})\pi_{ki}{\pi}_{km}$
$\displaystyle+\dfrac{1}{12}c^{2}_{k}\left\\{[L^{{}^{\prime\prime\prime}}(\bar{\bm{\uptheta}}^{(+)}_{k})+L^{{}^{\prime\prime\prime}}(\bar{\bm{\uptheta}}^{(-)}_{k})][\bm{\pi}_{k}\otimes\bm{\pi}_{k}\otimes\bm{\pi}_{k}])\;\pi_{km}\right\\}+\dfrac{\upepsilon^{(+)}_{k}-\upepsilon^{(-)}_{k}}{2c_{k}}\pi_{km}$

Based on the basic assumption of independence between $\upepsilon^{(\pm)}_{k}$
and $\bm{\pi}_{k}$ and the uncorrelation between $\pi_{ki}$ and $\pi_{kj}$,
$i\neq j$, we have:

$E\left[\dfrac{\upepsilon^{(+)}_{k}-\upepsilon^{(-)}_{k}}{2c_{k}}\bm{\pi}_{k}\right]=0;\;E(\pi^{2}_{km})=1;\;E(\pi_{ki}\pi_{km})=0.$

Then,

$\displaystyle E[\hat{g}_{km}(\hat{\bm{\uptheta}}_{k})]$
$\displaystyle=g_{m}(\hat{\bm{\uptheta}}_{k})E(\pi^{2}_{km})+\sum_{i=1,i\neq
m}^{p}g_{i}(\hat{\bm{\uptheta}}_{k})E(\pi_{ki}\pi_{km})$
$\displaystyle+\dfrac{c^{2}_{k}}{12}E\left\\{[L^{{}^{\prime\prime\prime}}(\bar{\bm{\uptheta}}^{(+)}_{k})+L^{{}^{\prime\prime\prime}}(\bar{\bm{\uptheta}}^{(-)}_{k})][\bm{\pi}_{k}\otimes\bm{\pi}_{k}\otimes\bm{\pi}_{k}]\;\pi_{km}\right\\}$
$\displaystyle=g_{m}(\hat{\bm{\uptheta}}_{k})+\dfrac{c^{2}_{k}}{12}E\left\\{[L^{{}^{\prime\prime\prime}}(\bar{\bm{\uptheta}}^{(+)}_{k})+L^{{}^{\prime\prime\prime}}(\bar{\bm{\uptheta}}^{(-)}_{k})][\bm{\pi}_{k}\otimes\bm{\pi}_{k}\otimes\bm{\pi}_{k}]\pi_{km}\right\\}$
$\displaystyle=g_{m}(\hat{\bm{\uptheta}}_{k})+\dfrac{1}{6}c^{2}_{k}E\left\\{[L^{{}^{\prime\prime\prime}}(\hat{\bm{\uptheta}}_{k})+O(c^{2}_{k})][\bm{\pi}_{k}\otimes\bm{\pi}_{k}\otimes\bm{\pi}_{k}]\;\pi_{km}\right\\}$
$\displaystyle=g_{m}(\hat{\bm{\uptheta}}_{k})+\dfrac{1}{6}c^{2}_{k}\left(\sum_{i,j,k}\dfrac{\partial^{3}L}{\partial\pi_{ki}\partial\pi_{kj}\partial\pi_{kl}}\;\pi_{ki}\pi_{kj}\pi_{kl}\right)\;\pi_{km}+O(c^{4}_{k})$
$\displaystyle=g_{m}(\hat{\bm{\uptheta}}_{k})+\dfrac{1}{6}c^{2}_{k}\left(\dfrac{\partial^{3}L(\hat{\bm{\uptheta}}_{k})}{\partial{\pi}_{km}^{3}}\;\pi^{4}_{km}+\sum_{i,j,l\neq
m}\dfrac{\partial^{3}L}{\partial\pi_{ki}\partial\pi_{kj}\partial\pi_{kl}}\;\pi_{ki}\pi_{kj}\pi_{kl}\pi_{km}\right)+O(c^{4}_{k})$

So, the $m$th bias term is in the form of:

$\displaystyle\begin{split}\bm{b}_{km}(\hat{\bm{\uptheta}}_{k})&=E[\hat{g}_{km}(\hat{\bm{\uptheta}}_{k})]-g_{m}(\hat{\bm{\uptheta}}_{k})\\\
&=\dfrac{1}{6}c^{2}_{k}\left(L^{{}^{\prime\prime\prime}}_{mmm}(\hat{\bm{\uptheta}}_{k})E(\pi^{4}_{km})+3\sum\limits_{i=1,i\neq
m}^{p}L^{{}^{\prime\prime\prime}}_{iim}(\hat{\bm{\uptheta}}_{k})E(\pi^{2}_{ki}\pi^{2}_{km})\right)+O(c^{4}_{k})\end{split}$

###  3.2 Convergence Conditions

This subsection presents conditions for convergence of the RDSA iteration:
$\bm{\uptheta}^{*}\rightarrow\bm{\uptheta}^{*}\;\text{a.s.}$. The following
conditions ensure convergence of $\hat{\bm{\uptheta}}_{k}$ to
$\bm{\uptheta}^{*}$ in the RDSA algorithm:

Consider all $k\geq K$ for some $K<\infty$. Suppose that for each $k$, there
exists some $\upalpha_{0},\upalpha_{1},\upalpha_{2}>0$, $\bm{\pi}_{k}$ is a
direction vector normalized so that $||\bm{\pi}_{k}||^{2}=p$ and symmetrically
distributed about 0 with $|\pi_{ki}^{4}|\leq\upalpha_{0}\;\text{a.s.
and}\;E|\pi_{ki}|^{2}\leq\upalpha_{1}.i=1,2,...,p$ For almost all
$\hat{\bm{\uptheta}}_{k}$, suppose that for $\forall\;\bm{\uptheta}$ in an
open neighborhood of $\hat{\bm{\uptheta}}_{k}$,
$L^{{}^{\prime\prime\prime}}(\bm{\uptheta})$ exists continuously with
individual elements satisfying
$|L^{{}^{\prime\prime\prime}}_{ijk}(\bm{\uptheta})|\leq\upalpha_{2}$. Similar
to the conditions for SPSA discussed in [15, Chap.7], if the derivatives of
$L$ are equicontinuous and bounded, then the RDSA algorithm will converge
almost surely under the following conditions:

  * •

A1: $a_{k},c_{k}>0,\;\forall k:a_{k}\rightarrow 0,c_{k}\rightarrow
0\;{\text{as}}\;k\rightarrow\infty,\;\sum\limits^{\infty}_{k=0}a_{k}=\infty,\;\sum\limits^{\infty}_{k=0}\left(\dfrac{a_{k}}{c_{k}}\right)^{2}<\infty.$

  * •

A2:
$\mathop{{\text{sup}}}\limits_{k}||\hat{\bm{\uptheta}}_{k}||<\infty\;\text{a.s.}$

  * •

A3: $\bm{\uptheta}^{*}$ is an asymptotically stable solution of the
differential equation: $d\bm{x}(t)/dt=-\bm{g}(\bm{x}).$

  * •

A4: Let
$D(\bm{\uptheta}^{*})=[x_{0}|\mathop{\text{lim}}\limits_{t\rightarrow\infty}\bm{x}(t|x_{0})=\bm{\uptheta}^{*}]$
where $\bm{x}(t|x_{0})$ denotes the solution to the differential equation of
$d\bm{x}(t)/dt=-\bm{g}(\bm{x})$ based on initial conditions $\bm{x}_{0}$.
There exists a compact $S\subseteq D(\bm{\uptheta}^{*})$ s.t
$\tilde{\bm{\uptheta}}\in S$ infinitely often for almost all sample points.

  * •

A5: $\forall
k,E|\upepsilon^{(\pm)}_{k}|^{2}\leq\upalpha_{0},E(\bm{\pi}_{k}\bm{\pi}^{T}_{k})=\bm{I}$,
and $E|\pi^{2}_{kl}L(\hat{\bm{\uptheta}}_{k}\pm
c_{k}{\bm{\pi}}_{k})^{2}|\leq\upalpha_{1},\;l=1,2,...,p.$

Proposition 1A: Let A1–A5 hold and for some $K<\infty$, each $k\geq K$,
${\\{\pi_{ki}\\}}$ is Gaussian (0,1) distributed i.i.d., then
$\hat{\bm{\uptheta}}_{k}\rightarrow\bm{\uptheta}^{*}$ as $k\rightarrow\infty$.

Proposition 1B: Let A1–A5 hold and for some $K<\infty$, each $k\geq K$,
${\\{\pi_{ki}\\}}$ are uniformly distributed on a $p$-dim sphere, then
$\hat{\bm{\uptheta}}_{k}\rightarrow\bm{\uptheta}^{*}$ as $k\rightarrow\infty$.

Proof: This proof applies for both 1A and 1B above. From the above subsection
in terms of bias in $\hat{\bm{g}}_{k}$ and the A1 condition, we know that:

$\displaystyle||b_{k}(\hat{\bm{\uptheta}}_{k})||<\infty\qquad\forall
k,\;b_{k}(\hat{\bm{\uptheta}}_{k})\rightarrow\bm{0}\;\;\text{a.s.}$ (4)

According to the inequality in [7, p. 315] , it holds that:

$\displaystyle P(\mathop{\text{sup}}\limits_{m\geq
k}||\sum\limits_{i=k}^{m}a_{i}\bm{e}_{i}||\geq\upeta)\leq\upeta^{-2}E||\sum\limits_{i=k}^{\infty}a_{i}\bm{e}_{i}||^{2}=\upeta^{-2}\sum\limits_{i=k}^{\infty}a^{2}_{i}E||\bm{e}_{i}||^{2}.$
(5)

By the definition of $\bm{e}_{k}(\hat{\bm{\uptheta}}_{k})$ and Condition A1,
A5, we would have:

$\mathop{\text{lim}}\limits_{k\rightarrow\infty}P(\mathop{\text{sup}}\limits_{m\geq
k}||\sum\limits_{i=k}^{m}a_{i}\bm{e}_{i}||\geq\upeta)=0\;\;\text{for\;any\;}\upeta>0.$

Then, by [5, Lemma 2.2.1], the above propositions hold.

###  3.3 Asymptotic Normality

This subsection mainly discusses the asymptotic normality for
$\hat{\bm{\uptheta}}_{k}$ in RDSA algorithm, especially, how the Fabian’s
Theorem in [5, Theorem 2.2] applies on RDSA algorithm based on its
generalization for SPSA shown in [12] and [6]. In order to satisfy the
conditions for Fabian’s Theorem, we strengthen the Condition A5 to:

A5′: For some $\updelta,\upalpha_{0},\upalpha_{1},\upalpha_{2}>0$ and $\forall
k,E|\upepsilon^{(\pm)}_{k}|^{2+\delta}\leq\upalpha_{0},E(\pi_{kl}L(\hat{\bm{\uptheta}}_{k}\pm
c_{k}\bm{\pi}_{k}))^{2+\updelta}\leq\upalpha_{1},E(\pi_{kl})^{4+\updelta}\leq\upalpha_{2},\;l=1,2,...,p,E(\bm{\pi}_{k}\bm{\pi}^{T}_{k})=\bm{I}$.

Then, the conditions for the asymptotic normality for
$\hat{\bm{\uptheta}}_{k}$ in RDSA can be stated as the following proposition:

Proposition 2: Assume that conditions A1–A4 and A5′ hold and $\upbeta>0$. Let
$\bm{H}(\cdot)$ denote the Hessian matrix for $L(\bm{\uptheta})$;
$\upsigma,\upphi$ be such that
$E(\upepsilon_{k}^{(+)}-\upepsilon_{k}^{(-)})^{2}\rightarrow\upsigma^{2}$ and
$E(\pi_{kl}^{4})\rightarrow\upphi,\;k\rightarrow\infty$. Suppose the gain
sequences are
$a_{k}=a/(k+1)^{\upalpha},c_{k}=c/(k+1)^{\upgamma},a>0,c>0,k=0,1,2,3,...$ and
set $0<\upalpha\leq 1,\upgamma\geq\upalpha/6,\upbeta=\upalpha-2\upgamma$. Let
$\bm{P}$ be an orthogonal matrix such that
$\bm{PH}(\bm{\uptheta}^{*})\bm{P}^{T}=a^{-1}\text{diag}(\uplambda_{1},...,\uplambda_{p})$.
Then

$\displaystyle
k^{\upbeta/2}(\hat{{\uptheta}}_{k}-{\uptheta}^{*})\stackrel{{\scriptstyle\text{dist}}}{{\longrightarrow}}N(\bm{\upmu},\bm{PMP}^{T}),\quad
k\rightarrow\infty,$ (6)

where
$\bm{M}=\dfrac{1}{4c^{2}}a^{2}\upsigma^{2}\text{diag}[1/(2\uplambda_{1}-\upbeta_{+}),....,1/(2\uplambda_{p}-\upbeta_{+})]$
with $\upbeta_{+}=\upbeta<2\mathop{\text{min}}\limits_{i}\uplambda_{i}$ if
$\upalpha=1$ and $\upbeta_{+}=0$ if $\upalpha<1$, and

$\displaystyle\bm{\upmu}=\left\\{\begin{array}[]{ccr}0&&{\text{if}\;3\upgamma-\upalpha/2>0}\\\
(a\bm{H}(\uptheta^{*})-\dfrac{1}{2}\upbeta_{+}\bm{I})^{-1}\bm{T}&&{\text{if}\;3\upgamma-\upalpha/2=0}\\\
\end{array}\right.,$ (9)

where the $l$th component of $\bm{T}$ is:

$\displaystyle\bm{T}^{\text{RD}}_{l}=-\dfrac{1}{6}ac^{2}\left[{\upphi}L^{{}^{\prime\prime\prime}}_{lll}(\bm{\uptheta}^{*})+3{\upsilon}\sum\limits_{i\neq
l}^{p}L^{{}^{\prime\prime\prime}}_{iil}(\bm{\uptheta}^{*})\right],$

and

$\displaystyle{\upphi}=E(\pi^{4}_{kl}),\;{\upsilon}=E(\pi^{2}_{kl}\pi^{2}_{km}),\;l\neq
m.$

Proof: The result will be shown if conditions (2.2.1), (2.2.2), and (2.2.3) of
[5, Theorem 2.2] hold. In the notation of [5, Theorem 2.2], we can rewrite the
formula (1) in the form of:

$\displaystyle\hat{\bm{\uptheta}}_{k+1}-\bm{\uptheta}^{*}=(\bm{I}-k^{-\upalpha}\bm{\Gamma}_{k})(\hat{\bm{\uptheta}}_{k}-\bm{\uptheta}^{*})+k^{-(\upalpha+\upbeta)/2}\bm{\Phi}_{k}\bm{V}_{k}+k^{-\upalpha-\upbeta/2}\bm{T}_{k},$

where

$\displaystyle\bm{\Gamma}_{k}$
$\displaystyle=k^{\upalpha}a_{k}\tilde{\bm{H}}_{k},\;\;\;\bm{\Phi}_{k}=-a\bm{I},$
$\displaystyle\bm{V}_{k}$
$\displaystyle=\dfrac{1}{k^{\upgamma}}[\hat{\bm{g}}_{k}(\hat{\bm{\uptheta}}_{k})-E(\hat{\bm{g}}_{k}(\hat{\bm{\uptheta}}_{k})],\;\;\bm{T}_{k}=-ak^{\upbeta/2}\bm{b}_{k}(\hat{\bm{\uptheta}}_{k}).$

Next, we are going to verify Fabian’s conditions.

As shown in the [4, p. 233], let
$\bm{\Gamma}_{k}=k^{\upalpha}a_{k}\tilde{\bm{H}}_{k}$ where the $i$th row of
$\tilde{\bm{H}}_{k}$ equals to the $i$th row of the Hessian matrix evaluated
at
$\bm{\uptheta}=(1-\uplambda_{i})\hat{\bm{\uptheta}}_{k}+\uplambda_{i}\bm{\uptheta}^{*}$
for some $\uplambda_{i}\in[0,1]$. Due to the continuity of
$\bm{H}(\bm{\uptheta})$ at $\bm{\uptheta}^{*}$ and
$k^{\upalpha}a_{k}\rightarrow a,k\rightarrow\infty$,
$a\bm{H}(\bm{\uptheta}^{*})$ would be a positive definite matrix, which
$\bm{\Gamma}_{k}$ would converge to.

Next, consider the convergence of $\bm{T}_{k}$:

$\displaystyle T_{kl}$
$\displaystyle=-\dfrac{1}{12}ac^{2}\dfrac{1}{k^{3\upgamma-\upalpha/2}}E\left\\{\pi_{kl}[L^{{}^{\prime\prime\prime}}(\bar{\bm{\uptheta}}^{(+)}_{k})+L^{{}^{\prime\prime\prime}}(\bar{\bm{\uptheta}}^{(-)}_{k})]\bm{\pi}_{k}\otimes\bm{\pi}_{k}\otimes\bm{\pi}_{k}\right\\}.$

If $3\upgamma-\upalpha/2>0$, we have ${T}_{kl}\rightarrow 0$ a.s.

If $3\upgamma-\upalpha/2=0$, use the fact that $L^{{}^{\prime\prime\prime}}$
is uniformly bounded near $\bm{\uptheta}^{*}$, we have:

$T_{kl}\rightarrow-\dfrac{ac^{2}}{6}\left\\{E(\pi^{4}_{kl})L^{{}^{\prime\prime\prime}}_{lll}(\bm{\uptheta}^{*})+3\sum\limits_{i\neq
l}^{p}L^{{}^{\prime\prime\prime}}_{iil}(\bm{\uptheta}^{*})E(\pi^{2}_{ki}\pi^{2}_{kl})\right\\}\;\text{a.s.}$

So, we can show that $\bm{T}_{k}$ converges a.s. for $3\upgamma-\upalpha/2\geq
0.$. Further, obviously,
$\bm{\Phi}_{k}=-a\bm{I}\rightarrow\bm{\Phi}=-a\bm{I}.$

Last, consider $\bm{V}_{k}$:

$\displaystyle\begin{split}&E(\bm{V}_{k}\bm{V}^{T}_{k})\\\
&=k^{-2\upgamma}E\bigg{\\{}\bm{\pi}_{k}(\bm{\pi}_{k})^{T}\left[\dfrac{L(\hat{\bm{\uptheta}}_{k}+c_{k}\bm{\pi}_{k})-L(\hat{\bm{\uptheta}}_{k}-c_{k}\bm{\pi}_{k})}{2ck^{-\upgamma}}\right]^{2}\\\
&+\bm{\pi}_{k}(\bm{\pi}_{k})^{T}\left[\dfrac{\upepsilon^{(+)}_{k}-\upepsilon^{(-)}_{k}}{2ck^{-\upgamma}}\right]\left[\dfrac{L(\hat{\bm{\uptheta}}_{k}+c_{k}\bm{\pi}_{k})-L(\hat{\bm{\uptheta}}_{k}-c_{k}\bm{\pi}_{k})}{2ck^{-\upgamma}}\right]\\\
&+\bm{\pi}_{k}(\bm{\pi}_{k})^{T}\left[\dfrac{\upepsilon^{(+)}_{k}-\upepsilon^{(-)}_{k}}{2ck^{-\upgamma}}\right]^{2}\bigg{\\}}\\\
&-k^{-2\upgamma}\left[\bm{g}(\hat{\bm{\uptheta}}_{k})+\bm{b}_{k}(\hat{\bm{\uptheta}}_{k})\right]\left[\bm{g}(\hat{\bm{\uptheta}}_{k})+\bm{b}_{k}(\hat{\bm{\uptheta}}_{k})\right]^{T}.\end{split}$
(10)

Similar to the proof (3.5) in [12], for sufficiently large $k$,
$L(\hat{\bm{\uptheta}}_{k}\pm c_{k}{\bm{\pi}_{k}})$ is uniformly bounded in
$\bm{\pi}_{k}$. Combined with condition A2 and Holder’s Inequality, it is
implied that the first and second terms in (10) will converge to 0 a.s. Also,
by conditions A1–A5, the fourth term will converge to 0 a.s. as well.

Then consider the third term:

$\displaystyle
E\left\\{\bm{\pi}_{k}(\bm{\pi}_{k})^{T}\left[\upepsilon^{(+)}_{k}-\upepsilon^{(-)}_{k}\right]^{2}\right\\}$
$\displaystyle=\int_{\Omega_{\bm{\pi}}}\bm{\pi}_{k}(\bm{\pi}_{k})^{T}E\left[\upepsilon^{(+)}_{k}-\upepsilon^{(-)}_{k}\right]^{2}dP_{\bm{\pi}}.$

where $\Omega_{\bm{\pi}}$ is the sample space generating the $\bm{\pi}_{k}$
and $P$ is the corresponding probability measure. Here we have the fact that
$E\left[\upepsilon^{(+)}_{k}-\upepsilon^{(-)}_{k}\right]^{2}\rightarrow\upsigma^{2}$
a.s. and $||\bm{\pi}_{k}||^{2}=p$. If we consider ${\\{\bm{\pi}_{ij}}\\}$ to
be independently generated with mean 0 and second moment 1 and then it leads
to the result that $E({\pi}_{kl}{\pi}_{km})=E({\pi}_{kl})E({\pi}_{km})=0$. If
we consider $\\{\pi_{k}\\}$ being generated from a $p$-dimensional sphere with
radius $\sqrt{p}$ uniformly, then for $l\neq m,\;{\pi}_{kl},{\pi}_{km}$ are
not independent, but uncorrelated. In this case, we still have
$E({\pi}_{kl}{\pi}_{km})=E({\pi}_{kl})E({\pi}_{km})=0$. Moreover,
$E({\pi}^{2}_{kl})=E({\pi}^{2}_{km})$, such that $E({\pi}^{2}_{kl})=1$. So, we
can have: $E[\bm{\pi}_{k}(\bm{\pi}_{k})^{T}]=\bm{I},\;\;\forall\;l\neq m$.
such that,

$E(\bm{V}_{k}\bm{V}^{T}_{k})\rightarrow\dfrac{1}{4}c^{-2}\upsigma^{2}\bm{I}.$

Now, the conditions (2.2.1) and (2.2.2) in [5, Theorem 2.2] have been shown.
Next, we are going to show the condition (2.2.3) also holds. By the definition
of $\bm{V}_{k}$, Markov Inequality, and the triangle inequality and the proof
shown in (3.8) in [12], for any $0<\upsilon<\updelta/2$, it holds that:

$\displaystyle||\bm{V}_{k}||^{2(1+\upsilon)}\leq$ $\displaystyle
2^{2(1+\upsilon)}k^{-2\upgamma(1+\upsilon)}[||\hat{\bm{g}}_{k}(\hat{\bm{\uptheta}}_{k})||^{2(1+\upsilon)}$
$\displaystyle+||\bm{b}_{k}(\hat{\bm{\uptheta}}_{k})||^{2(1+\upsilon)}+||\bm{g}_{k}(\hat{\bm{\uptheta}}_{k})||^{2(1+\upsilon)}].$

Since $\bm{g}(\hat{\bm{\uptheta}}_{k})$ and
$\bm{b}_{k}(\hat{\bm{\uptheta}}_{k})$ are uniformly bounded and
$L(\hat{\bm{\uptheta}}_{k}\pm c_{k}{\bm{\pi}}_{k})$ is uniformly bounded for
$\forall k>K$, we have

$\displaystyle
E[k^{-2\upgamma(1+\upsilon)}||\bm{b}_{k}(\hat{\bm{\uptheta}}_{k})||^{2(1+\upsilon)}]\rightarrow\bm{0},$
$\displaystyle
E[k^{-2\upgamma(1+\upsilon)}||\bm{g}_{k}(\hat{\bm{\uptheta}}_{k})||^{2(1+\upsilon)}]\rightarrow\bm{0}.$

as $k\rightarrow\infty$. Moreover, invoking A5′, Holder’s Inequality, similar
to the case in [12], it implies that
$E[||\hat{g}_{k}(\hat{\bm{\uptheta}}_{k})||^{2(1+\upsilon)}]=O(k^{2\upgamma(1+\upsilon)})$.
As a result, $E||\bm{V}_{k}||^{2(1+\upsilon)}=O(1)$, which shows that:

$\mathop{\rm{lim}}\limits_{k\rightarrow\infty}E\left(I_{||\bm{V}_{k}||^{2}\geq
rk^{\upalpha}}||\bm{V}_{k}||^{2}\right)=0,\qquad\forall r>0,$

where $I_{\\{.\\}}$ is the indicator function.

Now all the required conditions for Fabian’s Theorem in [5, Theorem 2.2] have
been verified. It is stated that $\hat{\bm{\uptheta}}_{k}$ in RDSA algorithm
is asymptotically normal.

Corollary 2A: Assume conditions in Proposition 2 hold and perturbation
$\\{{\bm{\pi}_{k}\\}}$ is Gaussian distributed $N(0,1)$ independently. Then
the $l$th component of $\bm{T}$ is:

$\bm{T}^{\text{RD}}_{l}=-\dfrac{1}{6}ac^{2}\left[3L^{{}^{\prime\prime\prime}}_{lll}(\bm{\uptheta}^{*})+3\sum\limits_{m\neq
l}^{p}L^{{}^{\prime\prime\prime}}_{llm}(\bm{\uptheta}^{*})\right].$

Corollary 2B: Assume conditions in Proposition 2 hold and perturbation
$\\{{\bm{\pi}_{k}}\\}$ is uniformly distributed on a $p$-dim sphere with
radius $\sqrt{p}$. Then the $l$th element of $\bm{T}$ is:

$\bm{T}^{\text{RD}}_{l}=-\dfrac{1}{6}ac^{2}\left[\dfrac{3p}{p+2}L^{{}^{\prime\prime\prime}}_{lll}(\bm{\uptheta}^{*})+\dfrac{3p}{p+2}\sum\limits_{m\neq
l}^{p}L^{{}^{\prime\prime\prime}}_{llm}(\bm{\uptheta}^{*})\right].$

###  3.4 Relative Accuracy between SPSA and RDSA

####  3.4.1 Introduction

This subsection discusses the MSE for SPSA and RDSA algorithm using the same
gain sequence. We mainly focus on several widely used distributions for
perturbations in SPSA and RDSA. We choose Bernoulli, and U-shape, especially
polynomial with high order for $\bm{\Delta}_{k}$. (Here we choose
$\Delta_{ki}\stackrel{{\scriptstyle\text{i.i.d}}}{{\sim}}x^{10},\,x\in[-1.17,1.17]$
to make it a valid probability density function) and choose Gaussian and
uniform-spherical for $\bm{\pi}_{k}$ in RDSA.

Notice that $\hat{{\bm{\uptheta}}}_{k}$ in the SPSA and RDSA algorithms have
the asymptotic normality in, respectively, [12] or formula (6) with different
values for the parameters. Here we give several notations to make the
parameters more concise:

$\displaystyle\begin{split}&\upphi=E[({\pi}^{4}_{ki})]\;\;;\;\,\qquad{\upsilon}=E({\pi}^{2}_{ki}{\pi}^{2}_{km}),\;i\neq
m\\\
&{\uprho}^{2}=E\left[({\Delta}_{ki})^{-2}\right];\;\;{\xi}^{2}=E\left[({\Delta}_{ki})^{2}\right].\end{split}$

Then we have:

$\displaystyle\bm{T}^{\text{RD}}_{l}$
$\displaystyle=-\dfrac{1}{6}ac^{2}\left[{\upphi}L^{{}^{\prime\prime\prime}}_{lll}(\bm{\uptheta}^{*})+3{\upsilon}\sum\limits_{m\neq
l}^{p}L^{{}^{\prime\prime\prime}}_{llm}(\bm{\uptheta}^{*})\right];$
$\displaystyle\bm{M}^{\text{RD}}$
$\displaystyle=\dfrac{1}{4c^{2}}a^{2}\upsigma^{2}\text{diag}\left(\dfrac{1}{2\uplambda_{1}-\beta_{+}},....,\dfrac{1}{2\uplambda_{p}-\beta_{+}}\right);$
$\displaystyle\bm{T}^{\text{SP}}_{l}$
$\displaystyle=-\dfrac{1}{6}ac^{2}{\xi}^{2}\left[L^{{}^{\prime\prime\prime}}_{lll}(\bm{\uptheta}^{*})+3\sum\limits_{m\neq
l}^{p}L^{{}^{\prime\prime\prime}}_{llm}(\bm{\uptheta}^{*})\right];$
$\displaystyle\bm{M}^{\text{SP}}$
$\displaystyle=\dfrac{1}{4c^{2}}a^{2}\upsigma^{2}{\uprho}^{2}\text{diag}\left(\dfrac{1}{2\uplambda_{1}-\upbeta_{+}},....,\dfrac{1}{2\uplambda_{p}-\beta_{+}}\right).$

Notice that tr$(\bm{M}_{\text{SP}})={\uprho}^{2}$tr$(\bm{M}_{\text{RD}}).$
Then, we can compute the asymptotic mean-square error by the definition:

$\displaystyle{\text{MSE}}=\bm{\upmu}^{T}\bm{\upmu}+{\text{tr}}(\bm{PMP}^{T})=\bm{\upmu}^{T}\bm{\upmu}+{\text{tr}}(\bm{M}).$

Asymptotically, the ratio of the MSE values of the two algorithms are in the
form of:

$\displaystyle\begin{split}\dfrac{\text{MSE}_{\text{RD}}}{\text{MSE}_{\text{SP}}}&\rightarrow\dfrac{\bm{\mu}^{T}_{\text{RD}}\bm{\mu}_{\text{RD}}+{\text{tr}}(\bm{PM}_{\text{RD}}\bm{P}^{T})}{\bm{\mu}^{T}_{\text{SP}}\bm{\mu}_{\text{SP}}+{\text{tr}}(\bm{PM}_{\text{SP}}\bm{P}^{T})}\\\
&={\left[{\upphi}\bm{u}_{1}+{\upsilon}\bm{u}_{2}\right]^{T}\bm{S}\left[{\upphi}\bm{u}_{1}+{\upsilon}\bm{u}_{2}\right]+{D}}{\left[\upxi^{2}\bm{u}_{1}+\upxi^{2}\bm{u}_{2}\right]^{T}\bm{S}\left[\upxi^{2}\bm{u}_{1}+\upxi^{2}\bm{u}_{2}\right]+{\uprho}^{2}{D}},\end{split}$
(11)

where

$\displaystyle\bm{u}_{1}=ac^{2}[L^{{}^{\prime\prime\prime}}_{111}(\bm{\uptheta}^{*}),...,L^{{}^{\prime\prime\prime}}_{ppp}(\bm{\uptheta}^{*})]^{T}/6;$
$\displaystyle\bm{u}_{2}=ac^{2}[3\sum\limits_{j\neq
1}L^{{}^{\prime\prime\prime}}_{jj1}(\bm{\uptheta}^{*}),...,3\sum\limits_{j\neq
p}L^{{}^{\prime\prime\prime}}_{jjp}(\bm{\uptheta}^{*})]^{T}/6;$
$\displaystyle\bm{S}=\left(a\bm{H}(\bm{\uptheta}^{*})-\dfrac{\upbeta_{+}}{2}\bm{I}\right)^{-2};$
$\displaystyle
D={\text{tr}}(\bm{M}_{\text{RD}})=\dfrac{a^{2}\upsigma^{2}}{4c^{2}}\sum\limits_{i=1}^{p}\dfrac{1}{2\uplambda_{i}-\upbeta_{+}}.$

Denote

$\displaystyle\begin{split}&{Q}_{1}=(\bm{u}_{1}+\bm{u}_{2})^{T}\bm{S}(\bm{u}_{1}+\bm{u}_{2}),\\\
&{Q}_{2}=(3\bm{u}_{1}+\bm{u}_{2})^{T}\bm{S}(3\bm{u}_{1}+\bm{u}_{2}).\end{split}$

The parameters ${\\{\upphi,\upsilon,\upxi^{2},\uprho^{2}}\\}$ and the MSE
values for several common used perturbations in SPSA and RDSA are listed as
follows:

Table 1: Key parameters for perturbations in RDSA and SPSA algorithms with different distributions Perturbation in RDSA | ${\upphi}$ | ${\upsilon}$
---|---|---
${\pi}_{ki}\stackrel{{\scriptstyle\text{i.i.d}}}{{\sim}}N(0,1)$ | 3 | 1
$\bm{\pi}_{k}\sim$ Uniformly Spherical | $\dfrac{3p}{p+2}$ | $\dfrac{p}{p+2}$
Perturbation in SPSA | ${\xi}^{2}$ | ${\rho}^{2}$
---|---|---
${\Delta}_{ki}\stackrel{{\scriptstyle\text{i.i.d}}}{{\sim}}$ Bernoulli ($\pm 1$) | 1 | 1
${\Delta}_{ki}\stackrel{{\scriptstyle\text{i.i.d}}}{{\sim}}x^{10}$, $x\in[-1.17,1.17]$ | 1.15 | 0.90
Table 2: MSE in SPSA and RDSA with different perturbations Algorithm | Distribution | MSE
---|---|---
SPSA | Bernoulli | ${Q}_{1}+D$
SPSA | U-shape ($x^{10}$) | $1.33{Q}_{1}+0.9D$
RDSA | Gaussian | ${Q}_{2}+D$
RDSA | Uniformly Spherical | $\dfrac{p}{p+2}{Q}_{2}+D$

From the two tables, we notice that when dimension $p$ is sufficiently large,
$\upphi,\upsilon$ of uniformly spherical are identical to those of the
Gaussian distribution. Also, as the order $d$ of $x^{d}$ in the U-shape
distribution increases, $\upxi^{2}$ of that distribution is close to that of
Bernoulli as well (the curve of the U-shape converges to two single points
$\pm 1$). Next, we are going to choose Bernoulli and Gaussian, respectively to
be the representative distribution of SPSA and RDSA and make comparisons
between them.

Proposition 3: If $2\bm{u}^{T}_{1}\bm{Su}_{1}+\bm{u}^{T}_{1}\bm{Su}_{2}\geq 0$
holds, then SPSA with Bernoulli distributed perturbation has a smaller MSE
value than RDSA with Gaussian perturbation.

Corollary 3A: If $|2\bm{u}_{1}|>|\bm{u}_{2}|$ holds component-wise and
$S_{ij}\geq 0$, then Proposition 3 holds.

Proof: $\bm{S}=(a\bm{H}(\bm{\uptheta}^{*})-\upbeta_{+}\bm{I}/2)^{-2}$ is a
positive-definite matrix by the fact that $\bm{H}(\bm{\uptheta}^{*})$ is
diagonalizable and
$\upbeta_{+}<2\text{min}\uplambda_{i}(a\bm{H}(\uptheta^{*}))$. By matrix
analysis, for any $\bm{u}_{1},\bm{u}_{2}\in\mathbb{R}^{p}$:

$|\bm{u}^{T}_{1}\bm{S}\bm{u}_{2}|\leq\sqrt{(\bm{u}^{T}_{1}\bm{Su}_{1})(\bm{u}^{T}_{2}\bm{Su}_{2})},$

So we would have,

$\displaystyle\begin{split}&2\bm{u}^{T}_{1}\bm{Su}_{1}+\bm{u}^{T}_{1}\bm{Su}_{2}\\\
&\geq
2\bm{u}^{T}_{1}\bm{Su}_{1}-\sqrt{(\bm{u}^{T}_{1}\bm{Su}_{1})(\bm{u}^{T}_{2}\bm{Su}_{2})}.\end{split}$
(12)

If $|2\bm{u}_{1}|\geq|\bm{u}_{2}|$ component-wise and $S_{ij}\geq 0$, then the
r.h.s of (12) would be larger than 0, so that
$2\bm{u}^{T}_{1}\bm{Su}_{1}+\bm{u}^{T}_{1}\bm{Su}_{2}\geq 0$, leading to,

$\displaystyle(3\bm{u}_{1}+\bm{u}_{2})^{T}\bm{S}(3\bm{u}_{1}+\bm{u}_{2})\geq(\bm{u}_{1}+\bm{u}_{2})^{T}\bm{S}(\bm{u}_{1}+\bm{u}_{2}).$

Consequently, the MSE of SPSA with Bernoulli-distributed perturbation will be
lower than that of RDSA with Gaussian-distributed perturbation.

Corollary 3B: If the loss function does not have cross-third derivative term,
i.e. $\bm{u}_{2}=\bm{0}$ along with $\upxi^{2}<\upphi,\uprho<1$, then SPSA has
a smaller MSE value than RDSA for any valid perturbation.

Proof: When $\bm{u}_{2}=\bm{0}$, we would have:

$\displaystyle\begin{split}&\text{MSE}_{\text{SPSA}}=\upxi^{4}\bm{u}^{T}_{1}\bm{Su}_{1}+\uprho^{2}D,\\\
&\text{MSE}_{\text{RDSA}}=\upphi^{2}\bm{u}^{T}_{1}\bm{Su}_{1}+D,\end{split}$
(13)

Then under the condition that $\upxi^{2}<\upphi,\uprho<1$, SPSA would have a
smaller MSE value.

####  3.4.2 Application on the skewed-quartic function

Consider the skewed quartic loss function:

$L(\bm{\uptheta})=\bm{\uptheta}^{T}\bm{B}^{T}\bm{B\uptheta}+0.1\sum\limits^{p}_{i=1}(\bm{B\uptheta})^{3}_{i}+0.01\sum\limits^{p}_{i=1}(\bm{B\uptheta})^{4}_{i}.$

where $\bm{\uptheta}=(t_{1},t_{2},...,t_{p})^{T}$. After calculation we obtain
that:

$\bm{u}_{1}=ac^{2}(L^{{}^{\prime\prime\prime}}_{111}(\bm{\uptheta}),...,L^{{}^{\prime\prime\prime}}_{ppp}(\bm{\uptheta}))^{T}/6|_{\bm{\uptheta}^{*}=\bm{0}}$=

$\displaystyle\dfrac{ac^{2}}{6}\left(\begin{array}[]{l}\dfrac{0.6}{p^{3}}+\dfrac{0.01}{p^{4}}\times
24\sum\limits^{p}_{i=1}t_{i}\\\ \dfrac{0.6}{p^{3}}\times
2+\dfrac{0.01}{p^{4}}\times
24\times(\sum\limits^{p}_{i=1}t_{i}+\sum\limits^{p}_{i=2}t_{i})\\\
\dfrac{0.6}{p^{3}}\times 3+\dfrac{0.01}{p^{4}}\times
24\times(\sum\limits^{p}_{i=1}t_{i}+\sum\limits^{p}_{i=2}t_{i}+\sum\limits^{p}_{i=3}t_{i})\\\
...\\\ \dfrac{0.6}{p^{3}}\times p+\dfrac{0.01}{p^{4}}\times
24\times(\sum\limits^{p}_{i=1}t_{i}+\sum\limits^{p}_{i=2}t_{i}+\sum\limits^{p}_{i=3}t_{i}+...+t_{p})\end{array}\right)_{\bm{\uptheta}^{*}=\bm{0}}=\dfrac{0.1\times
ac^{2}}{p^{3}}\left(\begin{array}[]{c}1\\\ 2\\\ 3\\\ ...\\\ p-1\\\
p\end{array}\right)$

$\bm{u}_{2}=ac^{2}(3\sum\limits^{p}_{j\neq
1}L^{{}^{\prime\prime\prime}}_{jj1}(\bm{\uptheta}),...,3\sum\limits^{p}_{j\neq
p}L^{{}^{\prime\prime\prime}}_{jjp}(\bm{\uptheta}))^{T}/6|_{\bm{\uptheta}^{*}=\bm{0}}$=

$\displaystyle\dfrac{ac^{2}}{2}\left(\begin{array}[]{l}(p-1)\left(\dfrac{0.6}{p^{3}}+\dfrac{0.01}{p^{4}}\times
24\sum\limits^{p}_{i=1}t_{i}\right)\\\
(p-1)\left(\dfrac{0.6}{p^{3}}+\dfrac{0.01}{p^{4}}\times
24\sum\limits^{p}_{i=1}t_{i}\right)+(p-2)\left(\dfrac{0.6}{p^{3}}+\dfrac{0.01}{p^{4}}\times
24\sum\limits^{p}_{i=2}t_{i}\right)\\\ ...\\\
(p-1)\left(\dfrac{0.6}{p^{3}}+\dfrac{0.01}{p^{4}}\times
24\sum\limits^{p}_{i=1}t_{i}\right)+...+\dfrac{0.6}{p^{3}}+\dfrac{0.01}{p^{4}}\times
24\sum\limits^{p}_{i=p-1}t_{i}\end{array}\right)$

$\displaystyle=\dfrac{0.3\times ac^{2}}{p^{3}}\left(\begin{array}[]{c}p-1\\\
(p-1)+(p-2)\\\ ...\\\ (p-1)+(p-2)+...+2+1\end{array}\right)$

$\bm{H}(\bm{\uptheta}^{*})=(L^{{}^{\prime\prime}}_{ij}(\bm{\uptheta}))_{ij}|_{\bm{\uptheta}^{*}=\bm{0}}$=

$\displaystyle\left(\begin{array}[]{ccccc}\dfrac{2}{p^{2}}&\dfrac{2}{p^{2}}&\dfrac{2}{p^{2}}&...&\dfrac{2}{p^{2}}\\\
\dfrac{2}{p^{2}}&\dfrac{2}{p^{2}}\times 2&\dfrac{2}{p^{2}}\times
2&...&\dfrac{2}{p^{2}}\times 2\\\ \dfrac{2}{p^{2}}&\dfrac{2}{p^{2}}\times
2&\dfrac{2}{p^{2}}\times 3&...&\dfrac{2}{p^{2}}\times 3\\\
\dfrac{2}{p^{2}}\times&\dfrac{2}{p^{2}}\times 2&\dfrac{2}{p^{2}}\times
3&...&\dfrac{2}{p^{2}}\times 4\\\ ...\\\
\dfrac{2}{p^{2}}&\dfrac{2}{p^{2}}\times 2&\dfrac{2}{p^{2}}\times
3&...&\dfrac{2}{p^{2}}\times
p\end{array}\right)=\dfrac{2}{p^{2}}\left(\begin{array}[]{cccccc}1&1&1&1&...&1\\\
1&2&2&2&...&2\\\ 1&2&3&3&...&3\\\ 1&2&3&4&...&4\\\ ...\\\
1&2&3&4&...&p\end{array}\right)$

In many practical cases, we set $\upalpha$ in the gain sequence
$a_{k}=a/(k+1)^{\upalpha}$ less than 1, so the $\upbeta^{+}=0$. Then

$\bm{S}=(a\bm{H}(\bm{\uptheta}^{*})-\dfrac{\upbeta^{+}}{2}\bm{I})^{-2}$=

$\displaystyle\dfrac{p^{4}}{4a^{2}}\left(\begin{array}[]{ccccccccc}5&-4&1&0&0&...&0&0&0\\\
-4&6&-4&1&0&...&0&0&0\\\ 1&-4&6&-4&1&...&0&0&0\\\ 0&1&-4&6&-4&...&0&0&0\\\
...\\\ 0&0&0&0&...&1&-4&6&-3\\\ 0&0&0&0&...&0&1&-3&2\end{array}\right)$

Then we could calculate the key terms:

$\displaystyle\begin{split}&\bm{u}^{T}_{1}\bm{Su}_{1}=\dfrac{0.01\times
c^{4}}{4a^{2}p^{2}};\\\ &\bm{u}^{T}_{1}\bm{Su}_{2}=0;\\\
&\bm{u}^{T}_{2}\bm{Su}_{2}=\dfrac{0.09\times
c^{4}(p-1)}{4a^{2}p^{2}}.\end{split}$

So, we will have:

$\displaystyle
2\bm{u}^{T}_{1}\bm{Su}_{1}+\bm{u}^{T}_{1}\bm{Su}_{2}=\dfrac{0.01\times
c^{4}}{2a^{2}p^{2}}.$ (14)

Since $p$ represents an integer larger than 1, the above term is always
positive. According to Proposition 3, it is known that the asymptotic MSE
value of SPSA with Bernoulli is smaller than the MSE for RDSA with Gaussian
(or other perturbation with fourth moment of its distribution being 3 and the
product of two separate second moments being 1). The corresponding numerical
experiment is shown in the next section.

Moreover, we can obtain the MSE value for each case:

$\displaystyle\begin{split}&Q_{1}=(\bm{u}_{1}+\bm{u}_{2})^{T}\bm{S}(\bm{u}_{1}+\bm{u}_{2})=\dfrac{0.09c^{4}p-0.08c^{4}}{4a^{2}p^{2}},\\\
&Q_{2}=(3\bm{u}_{1}+\bm{u}_{2})^{T}\bm{S}(3\bm{u}_{1}+\bm{u}_{2})=\dfrac{0.09c^{4}}{4a^{2}p},\\\
&D=\dfrac{\sigma^{2}}{8c^{2}}\sum\limits_{i}\dfrac{1}{\uplambda_{i}}\;\;\;\;(\uplambda_{i}\;\text{are\;the\;eigenvalues\;of\;matrix\;}a\bm{H}(\bm{\uptheta}^{*})).\end{split}$

From Table 2, we find that the asymptotic MSE value of Bernoulli will always
be smaller than that of Gaussian since $Q_{1}<Q_{2}$. Also, we find that the
MSE value of Uniform Spherical will be smaller than Gaussian with the
coefficient in front of $Q_{2}$ being less than 1. For the U-shape, there is
no certain conclusion for the relative MSE value since one contribution $D$ is
going down while the other one $Q_{1}$ is going up.

###  3.5 A More Straightforward Look

We can have a more straightforward understanding of the comparison of $Q_{1}$
and $Q_{2}$ by looking in a 2-dim way. Let $S,u_{1}$ and $u_{2}$ be scalars,
set $x=\sqrt{S}u_{1},y=\sqrt{S}u_{2}$, then consider $(3x+y)^{2}-(x+y)^{2}$:

$\displaystyle(3x+y)^{2}-(x+y)^{2}$ (15)
$\displaystyle=9x^{2}+6xy+y^{2}-x^{2}-2xy-y^{2}=8x^{2}+4xy.$ (16)

Under the following conditions, we will have (14) larger than 0.

$\left\\{\begin{array}[]{lr}y\geq-2x\;\;\;\text{if\;}x>0\\\
y<-2x\;\;\;\text{if\;}x<0\end{array}\right.$

The shadow region in the following figure shows the area satisfying the above
condition, which indicates a large probability, around 5/6 (more accurately,
1-2*arctan (0.5)) of $Q_{2}-Q_{1}>0$. In other words, there is a probability
around 5/6 that SPSA with Bernoulli outperforms RDSA with Gaussian.

Figure 1: Region of equation (14) being larger than 0

Then we are going to consider the high-dim case. Denote:

$\displaystyle\begin{split}&{Q}_{1}=(\bm{u}_{1}+\bm{u}_{2})^{T}\bm{S}(\bm{u}_{1}+\bm{u}_{2}),\\\
&{Q}_{2}=(3\bm{u}_{1}+\bm{u}_{2})^{T}\bm{S}(3\bm{u}_{1}+\bm{u}_{2}).\end{split}$

We would like to compare $Q_{2},Q_{1}$. Similar to the scalar case, set
$\bm{x}=\bm{S}^{1/2}\bm{u}_{1},\bm{y}=\bm{S}^{1/2}\bm{u}_{2}.$ Then consider:

$\displaystyle\bm{z}=Q_{2}-Q_{1}$ (17)
$\displaystyle=(3\bm{x}+\bm{y})^{T}(3\bm{x}+\bm{y})-(\bm{x}+\bm{y})^{T}(\bm{x}+\bm{y})$
(18) $\displaystyle=4(2\bm{x}^{T}\bm{x}+\bm{x}^{T}\bm{y})$ (19)
$\displaystyle=4\bm{x}^{T}(2\bm{x}+\bm{y})$ (20)

Here we consider a simple case to have a straightforward understanding. Set
$\bm{x},\bm{y}$ to be two independent random variables with uniform
symmetrical distribution (i.e.
$\bm{x}=(x_{1},x_{2},...,x_{p})^{T},\bm{y}=(y_{1},y_{2},...,y_{p})^{T},\;x_{i}\stackrel{{\scriptstyle\text{i.i.d}}}{{\sim}}U(-a,a),y_{i}\stackrel{{\scriptstyle\text{i.i.d}}}{{\sim}}U(-b,b)$)
and $\bm{S}$ to be the identical matrix. Then we have:

$\displaystyle
E(\bm{x}^{T}\bm{x})=E\left(\sum\limits^{p}_{i=1}x^{2}_{i}\right)=\sum\limits^{p}_{i=1}E(x^{2}_{i})=\sum\limits^{p}_{i=1}Var(x_{i})=\dfrac{a^{2}p}{3}.$
$\displaystyle
E(\bm{x}^{T}\bm{y})=E\left(\sum\limits^{p}_{i=1}x_{i}y_{i}\right)=\sum\limits^{p}_{i=1}E(x_{i})E(y_{i})=0.\;\;\;(\text{due\;to\;the\;independence})$
$\displaystyle
E(\bm{z})=8E(\bm{x}^{T}\bm{x})+4E(\bm{x}^{T}\bm{y})=\dfrac{8a^{2}p}{3}.$

Furthermore, we use one-sided version of Chebyshev inequality to see the
probability that $\bm{z}\leq 0$, which indicates $Q_{2}\leq Q_{1}$. The
inequality states that $P(X\geq c)\leq Var(X)/[Var(X)+c^{2}]$. In our case,
the inequality is in the form of:

$\displaystyle P[(-\bm{z}+E(\bm{z}))\geq E(\bm{z})]$
$\displaystyle\leq\dfrac{Var[-\bm{z}+E(\bm{z})]}{Var[-\bm{z}+E(\bm{z})]+[E(\bm{z})]^{2}}$
(21) $\displaystyle=\dfrac{Var(\bm{z})}{Var(\bm{z})+[E(\bm{z})]^{2}}.$ (22)

To simplify the calculation, we set
$x_{i}\stackrel{{\scriptstyle\text{i.i.d}}}{{\sim}}U(-a,a),y_{i}\stackrel{{\scriptstyle\text{i.i.d}}}{{\sim}}U(-a,a)$.
Then we have:

$\displaystyle Var(\bm{z})$
$\displaystyle=Var(8\bm{x}^{T}\bm{x})+Var(4\bm{x}^{T}\bm{y})+Cov(8\bm{x}^{T}\bm{x},4\bm{x}^{T}\bm{y})$
$\displaystyle=64Var\left(\sum\limits^{p}_{i=1}x^{2}_{i}\right)+16Var\left(\sum\limits^{p}_{i=1}x_{i}y_{i}\right)=\dfrac{64pa^{4}}{5}+\dfrac{16pa^{4}}{9}=16\times\dfrac{41pa^{4}}{45}.$

Plug it into the inequality (20), it looks like:

$\displaystyle P(\bm{z}\leq 0))$
$\displaystyle\leq\dfrac{Var(\bm{z})}{Var(\bm{z})+[E(\bm{z})]^{2}}$ (23)
$\displaystyle=\dfrac{16\times\dfrac{41pa^{4}}{45}}{16\times\dfrac{41pa^{4}}{45}+\dfrac{64p^{2}a^{4}}{9}}$
(24) $\displaystyle=\dfrac{41}{41+20p}.$ (25)

Inequality (23) shows that as the dimension $p$ increases, the probability of
$\bm{z}\leq 0$ decreases. Here we make a simulation by setting a range of
dimension from 1 to 10 and then take 100,000 independent trials to calculate
the probability, i.e. $P(\bm{z}\leq 0)$ with
$\bm{S}=\bm{I};\;\bm{u}_{1i}\stackrel{{\scriptstyle\text{i.i.d}}}{{\sim}}U(-100,100);\;\bm{u}_{2i}\stackrel{{\scriptstyle\text{i.i.d}}}{{\sim}}U(-100,100),\;i=1,2,...,p$.
The result is shown in the following Table 3, demonstrating a decreasing trend
when the dimension gets large. Also, it shows consistency with the conclusion
in the scalar case, namely, when $p=1,P(z\leq 0)\approx\dfrac{1}{6}$.

Table 3: Relationship Between Dimension and $P(\bm{z}\leq 0)$ Dim | 1 | 2 | 3 | 4 | 5
---|---|---|---|---|---
$\bm{P}(\bm{z}\leq 0)$ | 0.12546 | 0.0336 | 0.00932 | 0.00252 | 0.00071
Dim | 6 | 7 | 8 | 9 | 10
$\bm{P}(\bm{z}\leq 0)$ | 0.00024 | $9\times 10^{-5}$ | $3\times 10^{-5}$ | 0 | 0

As Table 3 indicates, if we consider those high-dimension cases, there will be
a tiny probability that $Q_{2}\leq Q_{1}$, i.e. SPSA with Bernoulli will
highly probably beat RDSA with Gaussian from the perspective of a smaller MSE.
However, it should be noticed that the above case is based on a simple
simulated setting, which aims at providing an intuitive view at the comparison
of Gaussian and Bernoulli.

In the next section, we will show how Proposition 3 works more specifically by
some numerical experiments.

##  4 Numerical Experiment

###  4.1 Numerical Example 1

Consider the loss function [3]:

$L(\uptheta)=||\uptheta||^{2}+\sum\limits_{i=1}^{p}e^{t_{i}/p},$

where $\bm{\uptheta}=[t_{1},t_{2},...,t_{p}]^{T}$ and the minimum occurs at
$\bm{\uptheta}^{*}$ with each component $-0.033$, such that
$L(\bm{\uptheta}^{*})=29.99994$. This is a convex function with third-
derivative being nonzero and cross third-derivative being zero. Namely, in our
case, we have $\bm{u}_{2}=\bm{0}$ and $\bm{u}_{1}>\bm{0}$. So according to the
Corollary in Proposition 3, SPSA with Bernoulli has lower MSE than RDSA with
Gaussian.

Let $p=30$ and take 100 independent trials, each with 3000 iterations. Results
are presented in Fig 2 and Table 4. We assume the noise term has distribution
$N(0,0.01)$ and gain sequence as
$a_{k}=0.05/(k+1)^{0.602},\,c_{k}=0.3/(k+1)^{0.101}$. Set the initial point
$\hat{\bm{\uptheta}}_{0}=[1,1,...,1]^{T}$. The results are shown below. The
lowest MSE is highlighted in bold in Table 4 and other numerical results
below.

Figure 2: Mean Square Error, computed as an average over 100 trials, plotted the last 200 iterations for different random perturbation distributions Table 4: MSE and confidence interval in SPSA and RDSA with different perturbations with same gain sequences in all algorithms Perturbation | MSE | 95% CI
---|---|---
Bernoulli SP | 0.01318 | [0.01281, 0.01354]
U-shape($\bm{x^{10}}$) SP | 0.01301 | [0.01266, 0.01336]
Gaussian RD | 0.01325 | [0.01290, 0.01359]
Uniformly Spherical RD | 0.01319 | [0.01285, 0.01355]

The experimental result shows consistency with the Proposition 3. The
confidence intervals show overlap among the four and then we conduct the two-
sample $t$-test between the Bernoulli and Gaussian cases and obtain the
$p$-value 0.35982, which is larger than the common threshold, 0.05, indicating
that the difference is not statistically significant. The reason is that in
this case, the $D$ value (shown in the Table 2) is much larger than the term
$\bm{u}_{1}^{T}\bm{Su}_{1}$ so that it dominates the whole MSE value. So the
MSE of SPSA with Bernoulli is slightly smaller than that of RDSA with
Gaussian. Moreover, we notice that the MSE of U-Shape distributed perturbation
is the smallest one which can also be explained by the 0.9 in front of the
value $D$.

###  4.2 Numerical Example 2

Consider the Ackley Function with dimension $p=30$:

$\displaystyle L(\bm{\uptheta})$
$\displaystyle=-a\times\text{exp}\left(-b\sqrt{\dfrac{1}{30}\sum\limits^{d}_{i=1}t^{2}_{i}}\right)-\text{exp}\left(\dfrac{1}{30}\sum\limits^{d}_{i=1}\text{cos}(ct_{i})\right)+a+e,$

where $a=20,b=0.2,c=2\pi,\,\bm{\uptheta}=[t_{1},t_{2},...,t_{30}]$ and
$L(\bm{\uptheta})=0$ at $\bm{\uptheta}^{*}=[1,1,...,1]^{T}.$ Here, we assume
the noise term has distribution $N(0,0.01)$ and gain sequence as
$a_{k}=0.02/(k+1+A)^{0.602},\;A=10,\,c_{k}=0.2/(k+1)^{0.101}$. Set the initial
point $\hat{\bm{\uptheta}}_{0}=[0.2,0.2,...,0.2]^{T}$.

In this case, there are no cross-third derivatives, i.e. $\bm{u}_{2}=\bm{0}$.
By Proposition 3, SPSA with Bernoulli would have smaller MSE value than RDSA
with Gaussian. Take 100 independent trials, each with 5000 iterations. The
results are shown in Fig 3 and Table 5:

Figure 3: Mean Square Error, computed as an average over 100 trials, plotted the last 500 iterations for different random perturbation distributions Table 5: MSE and its Confidence Interval in SPSA and RDSA with different perturbations with same gain sequences in all algorithms Perturbation | MSE | 95% CI
---|---|---
Bernoulli SP | 0.00041 | [0.00038, 0.00044]
U-shape ($x^{10}$) SP | 0.00043 | [0.00041, 0.00045]
Uniformly Spherical RD | 0.00054 | [0.00051, 0.00056]
Gaussian RD | 0.00057 | [0.00054, 0.00060]

The MSE value of SPSA with Bernoulli case is lower than that of RDSA with
Gaussian case, which is accordant with the Proposition 3. Notice that the
confidence interval of the two cases does not overlap and we take the two-
sample $t$-test between them and gain the $p$-value $5.96\times 10^{-13}$,
which indicates that the difference is significant.

###  4.3 Numerical Example 3

To give a more specific example. Consider the skewed-quartic function:

$L({\bm{\uptheta}})={\bm{\uptheta}}^{T}{\bm{B}}^{T}{\bm{B\uptheta}}+0.1\sum\limits^{p}_{i=1}({\bm{B\uptheta}})^{3}_{i}+0.01\sum\limits^{p}_{i=1}({\bm{B\uptheta}})^{4}_{i}.$

####  4.3.1 Optimal Gain Sequence Searching

We consider searching the ‘optimal gain sequence’ for $a_{k}$ and $c_{k}$ in
the gain sequence with the smallest MSE value, i.e. set
$a_{k}=a/(k+A+1)^{\alpha},\;c_{k}=c/(k+1)^{\gamma}$ and keep $A=10$. Then set
the search range of $a$ and $c$ to be 0.1 to 1, with interval 0.02 and pick
the sequence leading to the lowest MSE value. The ‘optimal gain sequence’ for
each case has been listed in Table 6. Particularly, in order to satisfy that
$\bm{\upmu}\neq\bm{0}$ in (9), we set $\upalpha=0.606,\upgamma=0.101$ in the
gain sequence.

Table 6: Optimal Gain Sequence for Algorithms with distinct Perturbations Perturbation | $a_{k}$ | $c_{k}$
---|---|---
Bernoulli SP | $0.12/(k+11)^{0.606}$ | $0.8/(k+1)^{0.101}$
U-shape ($x^{10}$) SP | $0.1/(k+11)^{0.606}$ | $0.48/(k+1)^{0.101}$
Uniformly Spherical RD | $0.1/(k+11)^{0.606}$ | $0.42/(k+1)^{0.101}$
Gaussian RD | 0.1/$(k+11)^{0.606}$ | $0.58/(k+1)^{0.101}$

####  4.3.2 Further Comparison

In order to conduct a further comparison of the performance of RDSA and SPSA
algorithm, as well as work in connection with our previous theoretical
calculation, we are going to apply the optimal gain sequence for each
distribution to the other three distributions and see the result. That is, in
1)-4) below, we, in turn, pick one of the pairs {$a_{k},c_{k}$} in Table 6 and
run all four methods at the chosen gain sequence. This will demonstrate each
method in its best configuration. Then we mainly focus on the Bernoulli and
Gaussian to see their comparison.

  * 1)

Apply the optimal gain sequence for Bernoulli, i.e. set
$a_{k}=0.12/((k+11)^{0.606}),c_{k}=0.8/((k+1)^{0.101}$ for the four methods.
The result is shown in Fig 4 and Table 7.

Figure 4: Mean Square Error, computed as an average over 100 trials each of
4000 iterations with the optimal gain for Bernoulli distributed perturbation
in SPSA, plotted the last 500 iterations for different random direction
distributions Table 7: MSE for Skewed Quartic Loss Function and its Confidence
Interval
(Gain Sequence $a_{k}=0.12/((k+11)^{0.606}),c_{k}=0.8/((k+1)^{0.101}$) Perturbation | MSE | 95% CI
---|---|---
Bernoulli SP | 6.3368 | [6.0959, 6.5778]
U-shape ($x^{10}$) SP | 6.3557 | [6.0536, 6.6579]
Uniformly Spherical RD | 6.6240 | [6.3064, 6.9416]
Gaussian RD | 6.6132 | [6.2637, 6.9626]

We find that the Bernoulli performs the best, then the U-shape, Uniformly
Spherical and Gaussian being the worst one. It is accordant with the previous
Table 2. Also, we find that the 95% confidence interval of Bernoulli and
Gaussian does not overlap. The two-sample test between Bernoulli and Gaussian
produces the $p$-value 0.0040573, indicating that the difference is
significant.

  * 2)

Apply the optimal gain sequence for U-shape, i.e. set
$a_{k}=0.1/((k+11)^{0.606}),c_{k}=0.48/((k+1)^{0.101}$ for the four methods.
The result is shown in Fig 5 and Table 8.

Figure 5: Mean Square Error, computed as an average over 100 trials each of
4000 iterations with the optimal gain for U-shape in SPSA, plotted the last
500 iterations for different random direction distributions Table 8: MSE for
Skewed Quartic Loss Function and its Confidence Interval
(Gain Sequence $a_{k}=0.1/((k+11)^{0.606}),c_{k}=0.48/((k+1)^{0.101}$) Perturbation | MSE | 95% CI
---|---|---
Bernoulli SP | 6.3784 | [6.1038, 6.6531]
U-shape($\bm{x^{10}}$) SP | 6.3540 | [6.0761, 6.6320]
Uniformly Spherical RD | 6.6257 | [6.3150, 6.9365]
Gaussian RD | 6.5604 | [6.2803, 6.8406]
  * 3)

Apply the optimal gain sequence for Uniform Spherical, i.e. set
$a_{k}=0.1/((k+11)^{0.606}),c_{k}=0.42/((k+1)^{0.101}$ for the four methods.
The result is shown in Fig 6 and Table 9.

Figure 6: Mean Square Error, computed as an average over 100 trials each of
4000 iterations with the optimal gain for Uniform Spherical in RDSA, plotted
the last 500 iterations for different random direction distributions Table 9:
MSE for Skewed Quartic Loss Function and its Confidence Interval
(Gain Sequence $a_{k}=0.1/((k+11)^{0.606}),c_{k}=0.42/((k+1)^{0.101}$) Perturbation | MSE | 95% CI
---|---|---
Bernoulli SP | 6.7323 | [6.4523,7.0123]
U-shape ($x^{10}$) SP | 6.7012 | [6.4212,6.9812]
Uniformly Spherical RD | 6.4436 | [6.1636,6.7236]
Gaussian RD | 7.0789 | [6.7889,7.3689]

We find that under such gain sequence, the Uniform Spherical performs the best
and Bernoulli still outperforms Gaussian, consistent with the Table 2. We then
conduct the two-sample $t$-test between Bernoulli and Gaussian and obtain the
$p$-value 0.121487, indicating that the difference is not statistically
significant.

  * 4)

Apply the optimal gain sequence for Gaussian, i.e. set
$a_{k}=0.1/((k+11)^{0.606}),c_{k}=0.58/((k+1)^{0.101}$ for the four methods.
The result is shown in Fig 7 and Table 10.

Figure 7: Mean Square Error, computed as an average over 100 trials each of
4000 iterations with the optimal gain for Gaussian distributed perturbation in
RDSA, plotted the last 500 iterations for different random direction
distributions Table 10: MSE for Skewed Quartic Loss Function and its
Confidence Interval
(Gain Sequence $a_{k}=0.1/((k+11)^{0.606}),c_{k}=0.58/((k+1)^{0.101}$) Perturbation | MSE | 95% CI
---|---|---
Bernoulli SP | 6.4928 | [6.2314, 6.7542]
U-shape ($x^{10}$) SP | 6.9915 | [6.7374, 7.2456]
Uniformly Spherical RD | 7.4873 | [7.1479, 7.8267]
Gaussian RD | 6.5534 | [6.2897, 6.8171]

In this case, although the optimal gain sequence of Gaussian has been applied
and the Gaussian reaches its lowest MSE, it is still smaller than the MSE of
Bernoulli, which can be a strong support for the theoretical conclusion shown
in Proposition 3. The two-sample test between Bernoulli and Gaussian produces
the $p$-value 0.2431523, indicating that the difference is not statistically
significant.

From above four experiments, all the four algorithms show their best
performance from the perspective of a lowest MSE when their optimal gain
sequence is applied. The tables show that no matter in which case, SPSA with
Bernoulli beats RDSA with Gaussian, showing consistency with Proposition 3.

##  5 Conclusion

This paper provides a strict derivation for several essential terms in the
RDSA algorithm and presents a comparison of the asymptotic MSE between RDSA
and SPSA. Furthermore, this paper presents specific conditions under which
SPSA with Bernoulli perturbation works better in the MSE sense than RDSA with
Gaussian perturbation, which may provide practitioners a general guidance in
the practical cases. In particular, it is found that SPSA outperforms RDSA
across the large majority of implementation cases with different gain
sequences and types of loss functions.

There are several directions for future research. First, we can consider extra
types of random directions with different distributions. Second, more work
related to robustness and relative efficiency would be useful. Third, since
this paper mainly focuses on Bernoulli and Gaussian cases, more general cases
could be considered for the comparison of MSE between SPSA and RDSA. It would
also be of interest to carry out a comparison such as that above for both
second-order or non-continuous versions of SPSA and RDSA, such as presented in
[16], [19], [17], [21], and [10], and constant-gain (step size) versions, such
as in [8], [20], and [11].

## References

  * [1] Shalabh Bhatnagar, Harshika Prasad, and L. A. Prashanth.  Stochastic recursive algorithms for optimization. Simultaneous perturbation methods, volume 434.  Jan 2013.
  * [2] A. Blakney and J. Zhu.  A comparison of the finite difference and simultaneous perturbation gradient estimation methods with noisy function evaluations.  In 2019 53rd Annual Conference on Information Sciences and Systems (CISS), pages 1–6, 2019.
  * [3] 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.
  * [4] Karla Hernández Cuevas.  Cyclic stochastic optimization: Generalizations, convergence, and applications in multi-agent systems, PhD Dissertation, Johns Hopkins University, https://arxiv.org/abs/1707.06700, 2017.
  * [5] Vaclav Fabian.  On asymptotic normality in stochastic approximation.  The Annals of Mathematical Statistics, 39(4):1327–1332, 1968.
  * [6] Karla Hernández and J. C. Spall.  Generalization of a result of Fabian on the asymptotic normality of stochastic approximation.  Automatica, 99:420 – 424, 2019.
  * [7] J.L.Doob.  Stochastic Process.  New York: Wiley, 1953.
  * [8] Itsuro Kajiwara, Keiichiro Furuya, and Shinichi Ishizuka.  Experimental verification of a real-time tuning method of a model-based controller by perturbations to its poles.  Mechanical Systems and Signal Processing, 107:396–408, 2018.
  * [9] H. J. Kusher and G. Yin.  Random directions methods in stochastic approximation.  In Proceedings of the 36th IEEE Conference on Decision and Control, pages 3430–3435 vol.4, 1997.
  * [10] L. A. Prashanth, S. Bhatnagar, M. Fu, and S. Marcus.  Adaptive system optimization using random directions stochastic approximation.  IEEE Transactions on Automatic Control, 62(5):2223–2238, 2017\.
  * [11] A. Ramaswamy.  DSPG: Decentralized simultaneous perturbations gradient descent scheme.  In 2020 28th Euromicro International Conference on Parallel, Distributed and Network-Based Processing (PDP), pages 54–62, Los Alamitos, CA, USA, Mar 2020. IEEE Computer Society.
  * [12] J. C. Spall.  Multivariate stochastic approximation using a simultaneous perturbation gradient approximation.  IEEE Transactions on Automatic Control, 37(3):332–341, 1992\.
  * [13] J. C. Spall.  Developments in stochastic optimization algorithms with gradient approximations based on function measurements.  In Proceedings of Winter Simulation Conference, pages 207–214, 1994\.
  * [14] J. C. Spall.  An overview of the simultaneous perturbation method for efficient optimization.  Johns Hopkins APL Technical Digest, 19:482–492, 1998.
  * [15] J. C. Spall.  Introduction to stochastic search and optimization: estimation, simulation, and control, volume 65.  John Wiley & Sons, 2005.
  * [16] J. C. Spall.  Feedback and weighting mechanisms for improving Jacobian estimates in the adaptive simultaneous perturbation algorithm.  IEEE Transactions on Automatic Control, 54(6):1216–1229, 2009.
  * [17] S. Sun and J. C. Spall.  SPSA method using diagonalized Hessian estimate.  In 2019 IEEE 58th Conference on Decision and Control (CDC), pages 4922–4927, 2019.
  * [18] J. Theiler and J. Alper.  On the choice of random directions for stochastic approximation algorithms.  IEEE Transactions on Automatic Control, 51(3):476–481, 2006\.
  * [19] L. Wang, J. Zhu, and J. C. Spall.  Mixed simultaneous perturbation stochastic approximation for gradient-free optimization with noisy measurements.  In 2018 Annual American Control Conference (ACC), pages 3774–3779, 2018.
  * [20] J. Zhu and J. C. Spall.  Stochastic approximation with nondecaying gain: Error bound and data-driven gain-tuning.  International Journal of Robust and Nonlinear Control, 30(15):5820–5870, 2020.
  * [21] J. Zhu, L. Wang, and J. C. Spall.  Efficient implementation of second-order stochastic approximation algorithms in high-dimensional problems.  IEEE Transactions on Neural Networks and Learning Systems, 31:3087–3099, 2020.