#  Theoretical Analysis of Explicit Averaging and Novel Sign Averaging in
Comparison-Based Search ††thanks: This study was partially supported by JSPS
KAKENHI under Grant Number 19H04179.

Daiki Morinaga D. Morinaga is with the Department of Computer Science,
University of Tsukuba, Tsukuba, Japan; RIKEN AIP, Japan (e-mail:
morinaga@bbo.cs.tsukuba.ac.jp). Youhei Akimoto  Y. Akimoto is with the
Institute of Systems and Information Engineering, University of Tsukuba,
Tsukuba, Japan; RIKEN AIP, Japan (e-mail: akimoto@cs.tsukuba.ac.jp).

###### Abstract

In black-box optimization, noise in the objective function is inevitable.
Noise disrupts the ranking of candidate solutions in comparison-based
optimization, possibly deteriorating the search performance compared with a
noiseless scenario. Explicit averaging takes the sample average of noisy
objective function values and is widely used as a simple and versatile noise-
handling technique. Although it is suitable for various applications, it is
ineffective if the mean is not finite. We theoretically reveal that explicit
averaging has a negative effect on the estimation of ground-truth rankings
when assuming stably distributed noise without a finite mean. Alternatively,
sign averaging is proposed as a simple but robust noise-handling technique. We
theoretically prove that the sign averaging estimates the order of the medians
of the noisy objective function values of a pair of points with arbitrarily
high probability as the number of samples increases. Its advantages over
explicit averaging and its robustness are also confirmed through numerical
experiments.

###### Index Terms:

comparison-based algorithm, explicit averaging, rank-based mechanism, sign
averaging, stochastic black-box optimization.

##  I Introduction

Evolutionary computation approaches in the continuous domain are intended to
approximate the optimal solution that minimizes objective function
$f:\mathbb{R}^{D}\to\mathbb{R}$ through queries $x\mapsto f(x)$ computed by
(possibly expensive) numerical simulations for each design variable $x$. With
the rapid development of computational performance, such approaches have been
increasingly exploited in practical and industrial applications, such as
geoscience [25], topology optimization [16], and machine learning [15]. A
major characteristic of most evolutionary computation approaches, such as the
evolution strategy with covariance matrix adaptation (CMA-ES) [20], and some
deterministic approaches, such as the Nelder–Mead method [26], is their
comparison-based principle. Instead of relying on the objective function
values, these approaches use the order of these values from the population of
design variables. Although the direct use of the objective function value has
been widely studied (e.g., in stochastic zeroth-order optimization [18]),
comparison-based approaches provide invariance against the strictly increasing
transformations of the objective function [8]. Under invariance, any
convergence statement for a single function can be generalized, both
theoretically and experimentally, to a class of functions that are
transformations of the original function. Furthermore, while the non-
differentiability or discontinuity of the objective function is often
troublesome for value-based approaches, it may be avoided if a strictly
increasing transformation exists for converting a troublesome problem into a
feasible one.

Real-world applications often entail the presence of noise in the $f$ values.
Consequently, each query $x\mapsto f(x;\bm{\epsilon})$ returns the objective
function value corrupted by noise $\bm{\epsilon}\in\mathrm{R}^{M}$, which is
independently sampled per query but remains concealed. Thus, the optimization
objective must consider the noise distribution, which is often defined as

$\operatornamewithlimits{argmin}_{x\in\mathbb{R}^{D}}\mathcal{A}[f(x;{\bm{\epsilon}})],$
(1)

where $\mathcal{A}$ denotes an optimization index. The choice of $\mathcal{A}$
depends on measures such as the expectation, percentile, or conditional value
at risk [17, 29]. In comparison-based search, the uncertainty of the objective
function affects the population rankings. Therefore, given a population of
design variables, $\\{x_{i}\\}_{i=1}^{\lambda}$, the search algorithms must
determine the population rankings with respect to the ground-truth objective
value, $\\{\mathcal{A}[f(x_{i};{\bm{\epsilon}})]\\}_{i=1}^{\lambda}$. However,
because search algorithms observe noisy objective function values, ground-
truth ranking cannot be achieved, thereby deteriorating the algorithm
performance.

Various noise-handling techniques have been proposed to address uncertainty in
evolutionary computation. These techniques can be classified into five groups
[29]: 1) estimating the ground truth by resampling and averaging the queries
(i.e., explicit averaging), 2) estimating the ground truth without using
explicit averaging, 3) enlarging the population size (implicit averaging), 4)
improving the search strategy, and 5) modifying the selection strategy.
Although no existing technique is predominant [7, 24, 6, 4, 28, 12], explicit
averaging is often employed in practice because it can be seamlessly
incorporated into most evolutionary computation approaches. Explicit averaging
resamples the noisy objective function values over $K\geqslant 1$ times per
$x_{i}$ and determines the population rankings, $\\{x_{i}\\}_{i=1}^{\lambda}$,
based on the average objective function values,
$\big{\\{}\frac{1}{K}\sum_{k=1}^{K}f(x_{i};\bm{\epsilon_{i,k}})\big{\\}}$.
Considering the strong law of large numbers, the averaged values converge to
$\\{\mathbb{E}[f(x_{i};{\bm{\epsilon}})]\\}_{i=1}^{\lambda}$ almost surely if
the distribution of the noisy objective function values has a finite mean.
Therefore, explicit averaging seems reasonable for handling noise in the
expectation as the optimization index.111 Explicit averaging has a main
advantage over implicit averaging. Assume that the ground-truth value
satisfies $f_{\mathrm{GT}}(x)=\mathbb{E}[f(x;{\bm{\epsilon}})]<\infty$.
Implicit averaging attempts to estimate the ranking of the objective function
values, $\\{f(x_{i};{\bm{\epsilon}})\\}_{i}^{\lambda}$, in population
$\\{x_{i}\\}_{i}^{\lambda}$. The algorithm estimation of the population
ranking can be evaluated using the Kendall rank correlation coefficient,
$\tau$, which measures the rank correlation between two sets of sample values
$\\{f(x_{i};{\bm{\epsilon}}_{i})\\}_{i}^{\lambda}$ and ground truth
$\\{\mathbb{E}[f(x_{i};{\bm{\epsilon}})]\\}_{i}^{\lambda}$. Coefficient $\tau$
is defined as
$\tau:=\frac{2}{\binom{\lambda}{2}}\cdot\sum_{i>j}^{\lambda}\mathbbm{1}{\left\\{(f(x_{i};{\bm{\epsilon}}_{i})-f(x_{j};{\bm{\epsilon}}_{j}))(\mathbb{E}[f(x_{i};{\bm{\epsilon}})]-\mathbb{E}[f(x_{j};{\bm{\epsilon}})])>0\right\\}}-1.$
Let
$p:=\Pr_{x,{\bm{\epsilon}}}\left[(f(x_{i};{\bm{\epsilon}}_{i})-f(x_{j};{\bm{\epsilon}}_{j})(\mathbb{E}[f(x_{i};{\bm{\epsilon}})]-\mathbb{E}[f(x_{j};{\bm{\epsilon}})])>0\right]<1$.
Considering the law of large numbers, $\lim_{\lambda\to\infty}\tau=2p-1<1$.
Therefore, implicit averaging cannot provide an exact estimate of ranking
$\tau=1$, even in the limit $\lambda\to\infty$. In contrast, with explicit
averaging,
$\frac{1}{K}\sum_{k}^{K}f(x;{\bm{\epsilon}}_{k})\to\mathbb{E}[f(x;{\bm{\epsilon}})]$
as $K\to\infty$ given the law of large numbers.

We theoretically investigate the effect of explicit averaging on comparison-
based algorithms. Existing theoretical studies on explicit averaging in
continuous black-box optimization often assume a finite variance for the noisy
objective function [10, 11, 9, 22, 28]. For finite variance, explicit
averaging reduces the standard deviation of the distribution of the noisy
objective function values by a factor of $1/\sqrt{K}$. However, it remains
unclear whether explicit averaging is effective for heavy-tailed distributions
of objective function values. In discrete optimization, Doerr and Sutton [14]
showed that averaging may fail when the noise follows a Cauchy distribution,
which does not have finite mean and variance. From these observations and the
strong law of large numbers holding even for infinite variance, we conjecture
that explicit averaging is effective if the noise distribution has a finite
mean, whereas it does not work otherwise. To the best of our knowledge, no
theoretical study has been conducted on the impact of explicit averaging on a
comparison-based search algorithm in the continuous domain for a heavy-tailed
noise distribution.

We consider two main research questions. 1) How effectively does averaging
with sample size $K$ work on comparison-based optimization algorithms when the
noise distribution does not have the $p$-th absolute moment as a finite value?
2) Can we design a noise-handling technique with theoretical grounds for a
noise distribution without a finite mean?

The main contributions of this study are summarized as follows:

  * •

In Section II, the effectiveness of explicit averaging with the stable noise
is addressed. Under 1 (linearity of $f(x;{\bm{\epsilon}})$ over
${\bm{\epsilon}}$) and 2 (stable noise), Theorem 3 reveals the range of stable
parameter $\alpha$ that determines whether explicit averaging works. For
$\alpha\in(1,2]$, explicit averaging is effective, while for $\alpha=1$, it
has no effect, and for $\alpha\in(0,1)$, explicit averaging is even harmful.
For stable noise, $\alpha$ determines the finiteness of the $p$-th absolute
moment.

  * •

In Section III, we propose a variant of explicit averaging called _sign
averaging_. In Theorem 9, we show that sign averaging is always effective
under different assumptions, namely, 3 (additivity of the median) and 4
(continuity of the distribution of $f(x;{\bm{\epsilon}})$). If the
distribution of $f(x;{\bm{\epsilon}})$ is symmetric, these two assumptions are
satisfied under 1 and 2.

##  II Averaging Analysis

We theoretically analyze the effect of increasing the sample size, $K$, by
averaging $\frac{1}{K}\sum_{k}^{K}f(x;{\bm{\epsilon}}_{k})$ on the _order
estimation probability_ (OEP), which is the probability of the estimate of the
order of two solutions, $x_{1}$ and $x_{2}$, agreeing with the ground truth.

###  II-A Stable Distribution

In this section, we assume that the noise follows a _stable distribution_ , as
we formally state below. The stable distribution is considered a
generalization of the well-studied distribution in noisy optimization, namely,
the _normal distribution_ , because the stable distribution is the only
possible distribution of the limit of the sum of a sequence of independent and
identically distributed (i.i.d.) random variables considering the generalized
central limit theorem [27].

Stable distribution ${\bm{S}}(\alpha,\beta,\gamma,\delta)$ is formally defined
as follows.

######  Definition 1 (Stable distribution
${\bm{S}}(\alpha,\beta,\gamma,\delta)$ in Definition 1.8 of [27]).

Let $\alpha\in(0,2]$, $\beta\in[-1,1]$, $\gamma\in(0,\infty)$, and
$\delta\in\mathbb{R}$. In addition, let random variable $X$ with support
$\mathbb{R}$ have the following characteristic function:

$\mathbb{E}[\exp(itX)]=\\\
\begin{cases}\exp\left(-\gamma^{\alpha}|t|^{\alpha}\left(1-\mathrm{i}\mkern
1.0mu\beta\tan\left(\frac{\pi\alpha}{2}\right)\mathrm{sign}(t)\right)+\mathrm{i}\mkern
1.0mu\delta t\right)&\text{ for }\alpha\neq 1;\\\
\exp\left(-\gamma|t|\left(1+\mathrm{i}\mkern
1.0mu\beta\frac{2}{\pi}\mathrm{sign}(t)\log|t|\right)+\mathrm{i}\mkern
1.0mu\delta t\right)&\text{ for }\alpha=1.\end{cases}$ (2)

The probability distribution of $X$ is called a stable distribution and
denoted by ${\bm{S}}(\alpha,\beta,\gamma,\delta)$, where $\mathrm{i}\mkern
1.0mu$ is the imaginary unit.

Parameters $\alpha$, $\beta$, $\gamma$, and $\delta$ are often referred to as
stability parameter, skewness parameter, scale parameter, and location
parameter, respectively.

Stable distributions include well-known distributions. For example, ${\bf
S}(\alpha=2,\beta,\gamma,\delta)$ is a normal distribution with mean $\delta$
and variance $2\gamma^{2}$, while ${\bf S}(\alpha=1,\beta=0,\gamma,\delta)$ is
a Cauchy distribution with mode $\delta$, and the Levy distribution is
obtained as ${\bf S}(\alpha=1/2,\beta=1,\gamma,\delta)$. The skewness
parameter controls the symmetricity of the distribution. If $\beta=0$, the
distribution is symmetric. Note that if $\alpha=2$, $\beta$ is ignored and
${\bf S}(2,\beta,\gamma,\delta)$ is always symmetric.

An important property of this distribution is integrability. The variance
exists as a finite value if and only if $\alpha=2$, that is, a normal
distribution. For $0<\alpha<2$, the $p$-th absolute moment,
$\mathbb{E}[|X|^{p}]$, is finite only for $0<p<\alpha$ [27]. Therefore, not
even the mean is defined for $0<\alpha\leqslant 1$.

The stable distribution is said to be _stable_ because a linear combination of
two random variables following stable distributions with the same $\alpha$ is
again a random variable following a stable distribution with the same
$\alpha$. This property is formally stated as follows.

######  Proposition 1 (Linear transformation of stable distribution in
Proposition 1.17 of [27]).

The following properties hold for stable distributions:

  1. 1.

If $X\sim{\bm{S}}(\alpha,\beta,\gamma,\delta)$, then for any
$a,b\in\mathbb{R}$,

$aX+b\sim\\\
\begin{cases}{\bm{S}}(\alpha,\mathrm{sign}(a)\cdot\beta,|a|\gamma,a\delta+b)&\alpha\neq
1;\\\
{\bm{S}}\left(1,\mathrm{sign}(a)\cdot\beta,|a|\gamma,a\delta+b-\frac{2}{\pi}\beta\gamma
a\log|a|\right)&\alpha=1.\end{cases}$ (3)

  2. 2.

If $X_{1}\sim{\bm{S}}(\alpha,\beta_{1},\gamma_{1},\delta_{1})$ and
$X_{2}\sim{\bm{S}}(\alpha,\beta_{2},\gamma_{2},\delta_{2})$ are independent,
then $X_{1}+X_{2}\sim{\bm{S}}(\alpha,\beta,\gamma,\delta)$, where

$\beta=\frac{\beta_{1}\gamma_{1}^{\alpha}+\beta_{2}\gamma_{2}^{\alpha}}{\gamma_{1}^{\alpha}+\gamma_{2}^{\alpha}},\enspace\gamma=\left(\gamma_{1}^{\alpha}+\gamma_{2}^{\alpha}\right)^{\frac{1}{\alpha}},\enspace\delta=\delta_{1}+\delta_{2}.$
(4)

###  II-B Problem Statement

We assume that $f(x;{\bm{\epsilon}})$ is linear with respect to
${\bm{\epsilon}}$ and that each element of
${\bm{\epsilon}}=[\epsilon_{1},\dots,\epsilon_{M}]^{\mathrm{T}}$ follows a
stable distribution, as formally stated below.

######  Assumption 1.

Function $f:\mathbb{R}^{D}\times\mathbb{R}^{M}\to\mathbb{R}$ is defined as

$f(x;{\bm{\epsilon}})=h(x)+\sum_{m}^{M}g_{m}(x)\cdot\epsilon_{m}$ (5)

for some $h:\mathbb{R}^{D}\to\mathbb{R}$ and
$g_{m}:\mathbb{R}^{D}\to\mathbb{R}$ with $m=1,2,\ldots,M$.

######  Assumption 2.

Noise vector ${\bm{\epsilon}}=[\epsilon_{1},\dots,\epsilon_{M}]^{\mathrm{T}}$
consists of random variables that follow a stable distribution with a common
$\alpha\in(0,2]$. The distribution of the $m$–th element, $\epsilon_{m}$, of
${\bm{\epsilon}}$ is ${\bf S}(\alpha,\beta_{m},\gamma_{m},\delta_{m})$.

The problem classes defined by 1 and 2 include frequently studied problem
instances. A function with additive noise
$f(x;{\bm{\epsilon}})=h(x)+\epsilon_{1}$ is included when $M=1$ and
$g_{1}(x)=1$. A function with multiplicative noise
$f(x;{\bm{\epsilon}})=h(x)(1+\epsilon_{1})$ is included when $M=1$ and
$g_{1}(x)=h(x)$. Noise is often assumed to follow a normal distribution, which
is an instance of a stable distribution.

The studied optimization index is given by

$\mathcal{A}[f(x;{\bm{\epsilon}})]=f(x;\Delta),$ (6)

where $\Delta=[\delta_{1},\dots,\delta_{M}]^{\mathrm{T}}$ collects the
location parameter, $\delta_{m}$, of $\epsilon_{m}$. This choice is explained
as follows. Arguably, the most common optimization index is expected value
$\mathbb{E}_{\bm{\epsilon}}[f(x;{\bm{\epsilon}})]$, and explicit averaging
likely minimizes the expected value. However, this is undefined for
$\alpha\leqslant 1$ unless $g_{m}(x)=0$ for all $m=1,\dots,M$. If
$\alpha\in(1,2]$, we have
$f(x;\Delta)=h(x)+\sum_{m}^{M}g_{m}(x)\cdot\delta_{m}=\mathbb{E}_{{\bm{\epsilon}}}[f(x;{\bm{\epsilon}})]$
[27]. Optimization index $f(x;\Delta)$ is well-defined for $\alpha\in(0,1]$,
whereas $\mathbb{E}_{\bm{\epsilon}}[f(x;{\bm{\epsilon}})]$ is undefined.
Additionally, if $\beta_{m}=0$ for all $m=1,\dots,M$, $f(x;\Delta)$
corresponds to the median of $f(x;\Delta)$ for all $\alpha\in(0,2]$.
Therefore, our objective can be also understood as minimizing the median of
the objective function value when $\beta_{m}=0$ for all $m=1,\dots,M$.
Hereafter,
$\nabla_{\bm{\epsilon}}f(x;{\bm{\epsilon}})=[g_{1}(x),\ldots,g_{M}(x)]^{\mathrm{T}}$
is denoted by $\nabla_{\bm{\epsilon}}f(x)$ to avoid confusion because it is
independent of ${\bm{\epsilon}}$.

###  II-C OEP

Comparison-based search algorithms determine the _order_ of any pair of
candidate solutions $(x_{1},x_{2})$. For a stochastic objective function, the
ground-truth order determined by $f(x_{1};\Delta)$ and $f(x_{2};\Delta)$ is
unknown and needs to be estimated from finite samples
$\\{f(x_{1};{\bm{\epsilon}}_{1,k})\\}_{k=1}^{K}$ and
$\\{f(x_{2};{\bm{\epsilon}}_{2,k})\\}_{k=1}^{K}$, where
${\bm{\epsilon}}_{i,k}$ are all i.i.d. copies of ${\bm{\epsilon}}$ satisfying
2. A common estimation approach is _explicit averaging_ , in which the order
of $f(x_{1};\Delta)$ and $f(x_{2};\Delta)$ is estimated by using the order of
$\frac{1}{K}\sum_{k=1}^{K}f(x_{1};{\bm{\epsilon}}_{1,k})$ and
$\frac{1}{K}\sum_{k=1}^{K}f(x_{2};{\bm{\epsilon}}_{1,k})$ with $K$ possibly
adapted during search.

We analyze the OEP, which is the probability of correctly estimating the
orders of $f(x_{1};\Delta)$ and $f(x_{2};\Delta)$. The OEP is formally defined
as follows. Let $\mathrm{sign}(\cdot)$ denote the sign function whose output
($-1,0$, or $1$) is the sign of its argument.

######  Definition 2 (OEP for $(x_{1},\Delta)$ and $(x_{2},\Delta)$).

Let $\widehat{\eta}_{{\bm{\epsilon}}}(x_{1},x_{2})\in\\{-1,0,1\\}$ be an
estimate of

$\eta_{\Delta}(x_{1},x_{2}):=\mathrm{sign}(f(x_{1};\Delta)-f(x_{2};\Delta)).$
(7)

The _OEP_ for $(x_{1},\Delta)$ and $(x_{2},\Delta)$ is given by

$\Pr_{\bm{\epsilon}}[\widehat{\eta}_{\bm{\epsilon}}(x_{1},x_{2})=\eta_{\Delta}(x_{1},x_{2})].$
(8)

The order estimate using explicit averaging is defined as

$\widehat{\eta}_{\epsilon}^{\mathrm{AVE}}(x_{1},x_{2}):=\mathrm{sign}\left(\frac{1}{K}\sum_{k}^{K}f(x_{1};{\bm{\epsilon}}_{1,k})-\frac{1}{K}\sum_{k}^{K}f(x_{2};{\bm{\epsilon}}_{2,k})\right).$
(9)

Under 1 and 2, the argument of $\mathrm{sign}(\cdot)$, namely,
$\frac{1}{K}\sum_{k=1}^{K}f(x_{1};{\bm{\epsilon}}_{1,k})-\frac{1}{K}\sum_{k=1}^{K}f(x_{2};{\bm{\epsilon}}_{2,k})$,
is stably distributed as stated below. Let the diagonal matrix whose
$(n,n)$-th element is $X_{n}$ be denoted by
$\mathrm{diag}(X_{1},\ldots,X_{N})$ and $\lVert
X\rVert_{\alpha}=(\sum_{n}^{N}|X_{n}|^{\alpha})^{1/\alpha}$ for
$X=[X_{1},\ldots,X_{N}]^{\mathrm{T}}$ be the $\alpha$-norm for
$\alpha\geqslant 1$. Let the indicator function, which is 1 if proposition $P$
is true and 0 otherwise, be denoted as $\mathbbm{1}{\left\\{P\right\\}}$. The
proof is presented in Section A-A.

######  Lemma 2 (Distribution of the difference of two averages).

Let $f(x;{\bm{\epsilon}})$ satisfy 1 and
${\bm{\epsilon}}=[\epsilon_{1},\ldots,\epsilon_{M}]^{\mathrm{T}}$ satisfy 2.
In addition, let ${\bm{\epsilon}}_{i,k}(i=1,2$ and $k=1,\ldots,K)$ be i.i.d.
copies of ${\bm{\epsilon}}$,
$\Gamma=\mathrm{diag}(\gamma_{1},\ldots,\gamma_{M})$, and
$\Delta=[\delta_{1},\ldots,\delta_{M}]^{\mathrm{T}}$. Define

$\displaystyle\beta^{\prime}(x)$
$\displaystyle=\frac{\sum_{m}^{M}\beta_{m}\cdot\mathrm{sign}(g_{m}(x))\cdot|g_{m}(x)\gamma_{m}|^{\alpha}}{\gamma^{\prime}(x)^{\alpha}},$
(10) $\displaystyle\gamma^{\prime}(x)$
$\displaystyle=\lVert\Gamma\nabla_{\bm{\epsilon}}f(x)\rVert_{\alpha},$ (11)
$\displaystyle\delta^{\prime}(x)$
$\displaystyle=-\frac{2}{\pi}\cdot\sum_{m}^{M}\beta_{m}\gamma_{m}g_{m}(x)\cdot\log|g_{m}(x)|,$
(12)

and

$\displaystyle\beta^{\prime\prime}(x_{1},x_{2})$
$\displaystyle=\frac{\beta^{\prime}(x_{1})\gamma^{\prime}(x_{1})^{\alpha}-\beta^{\prime}(x_{2})\gamma^{\prime}(x_{2})^{\alpha}}{\gamma^{\prime}(x_{1})^{\alpha}+\gamma^{\prime}(x_{2})^{\alpha}},$
(13) $\displaystyle\gamma^{\prime\prime}(x_{1},x_{2})$
$\displaystyle=(\gamma^{\prime}(x_{1})^{\alpha}+\gamma^{\prime}(x_{2})^{\alpha})^{\frac{1}{\alpha}},$
(14) $\displaystyle\delta^{\prime\prime}(x_{1},x_{2})$
$\displaystyle=\delta^{\prime}(x_{1})-\delta^{\prime}(x_{2}).$ (15)

It follows that

$\frac{1}{K}\sum_{k}^{K}f(x_{1};{\bm{\epsilon}}_{1,k})-\frac{1}{K}\sum_{k}^{K}f(x_{2};{\bm{\epsilon}}_{2,k})\\\
\sim{\bm{S}}\Big{(}\alpha,\beta^{\prime\prime}(x_{1},x_{2}),K^{\frac{1}{\alpha}-1}\cdot\gamma^{\prime\prime}(x_{1},x_{2}),\\\
f(x_{1};\Delta)-f(x_{2};\Delta)+\delta^{\prime\prime}(x_{1},x_{2})\cdot\mathbbm{1}{\left\\{\alpha=1\right\\}}\Big{)}.$
(16)

The following main result allows to evaluate the OEP for
$\widehat{\eta}_{\epsilon}^{\mathrm{AVE}}(x_{1},x_{2})$. The proof is
presented in Section A-B.

######  Theorem 3 (OEP over stable distribution).

Let $f(x;{\bm{\epsilon}})$ satisfy 1 and
${\bm{\epsilon}}=[\epsilon_{1},\ldots,\epsilon_{M}]^{\mathrm{T}}$ satisfy 2.
In addition, let ${\bm{\epsilon}}_{i,k}(i=1,2$ and $k=1,\ldots,K)$ be i.i.d.
copies of ${\bm{\epsilon}}$,
$\Gamma=\mathrm{diag}(\gamma_{1},\ldots,\gamma_{M})$, and
$\Delta=[\delta_{1},\ldots,\delta_{M}]^{\mathrm{T}}$. Define

$\epsilon^{\mathrm{AVE}}\sim{\bm{S}}\Bigg{(}\alpha,\beta^{\prime\prime}(x_{1},x_{2}),1,\\\
\left(\frac{\delta^{\prime\prime}(x_{1},x_{2})}{\gamma^{\prime\prime}(x_{1},x_{2})}+\frac{2}{\pi}\beta^{\prime\prime}(x_{1},x_{2})\log(\gamma^{\prime\prime}(x_{1},x_{2}))\right)\cdot\mathbbm{1}{\left\\{\alpha=1\right\\}}\Bigg{)}.$
(17)

For any pair $(x_{1},x_{2})$ with $\eta_{\Delta}(x_{1},x_{2})\neq 0$, we have

$\Pr_{\bm{\epsilon}}\left[\widehat{\eta}_{\epsilon}^{\mathrm{AVE}}(x_{1},x_{2})=\eta_{\Delta}(x_{1},x_{2})\right]\\\
=\Pr_{\bm{\epsilon}}\left[-\eta_{\Delta}(x_{1},x_{2})\cdot\epsilon^{\mathrm{AVE}}<\frac{K^{1-\frac{1}{\alpha}}\cdot|f(x_{1};\Delta)-f(x_{2};\Delta)|}{\gamma^{\prime\prime}(x_{1},x_{2})}\right].$
(18)

If $\eta_{\Delta}(x_{1},x_{2})=0$, then
$\Pr_{\bm{\epsilon}}\left[\widehat{\eta}_{\epsilon}^{\mathrm{AVE}}(x_{1},x_{2})=\eta_{\Delta}(x_{1},x_{2})\right]=0$.

The implications of Theorem 3 are as follows, noting that the distribution of
$\epsilon^{\mathrm{AVE}}$ is independent of $K$.

  * •

For $1<\alpha\leqslant 2$, the OEP approaches $1$ as $K$ increases because the
upper bound in the probability increases to infinity. However, the increase
rate of the OEP with $K$ is slower for a smaller $\alpha$.

  * •

For $\alpha=1$, explicit averaging has no effect on the OEP. A high $K$ simply
requires more $f$ queries.

  * •

For $0<\alpha<1$, the OEP decreases as $K$ increases because the upper bound
in the probability approaches 0.

The success of explicit averaging for $1<\alpha\leqslant 2$ can be understood
considering the law of large numbers. Average
$\frac{1}{K}\sum_{k}^{K}f(x;{\bm{\epsilon}}_{k})$ converges to expected value
$\mathbb{E}_{{\bm{\epsilon}}}[f(x;{\bm{\epsilon}})]=f(x;\Delta)$ if it exists
as a finite value. If the expected value does not exist ($0<\alpha\leqslant
1$), explicit averaging has a negative impact on the OEP, unless $\alpha=1$,
for which no impact occurs other than the increased cost of more $f$ queries
for a greater $K$.

To further clarify the OEP, we assume that the distribution is symmetric, that
is, $\beta_{m}=0$ for all $m=1,\dots,M$. In this case, both
$\beta^{\prime\prime}(x_{1},x_{2})$ and $\delta^{\prime\prime}(x_{1},x_{2})$
in Lemma 2 are zero for any pair $(x_{1},x_{2})$, and the expression of the
OEP is simplified. The proof is presented in Section A-C.

######  Corollary 4 (OEP over ${\bf S}(\alpha,0,\gamma,\delta)$).

Assume $\beta_{m}=0$ for all $m=1,\ldots,M$. Let $\phi_{\alpha}(\cdot)$ be the
cumulative distribution function of ${\bf S}(\alpha,0,1,0)$. For any pair
$(x_{1},x_{2})$ with $\eta_{\Delta}(x_{1},x_{2})\neq 0$, we have

$\Pr_{\bm{\epsilon}}\left[\widehat{\eta}_{\epsilon}^{\mathrm{AVE}}(x_{1},x_{2})=\eta_{\Delta}(x_{1},x_{2})\right]\\\
=\phi_{\alpha}\left(\frac{K^{1-\frac{1}{\alpha}}|f(x_{1};\Delta)-f(x_{2};\Delta)|}{\left(\|\Gamma\nabla_{\bm{\epsilon}}f(x_{1})\|_{\alpha}^{\alpha}+\|\Gamma\nabla_{\bm{\epsilon}}f(x_{2})\|_{\alpha}^{\alpha}\right)^{\frac{1}{\alpha}}}\right).$
(19)

If $\eta_{\Delta}(x_{1},x_{2})\neq 0$, then
$\Pr_{\bm{\epsilon}}\left[\widehat{\eta}_{\epsilon}^{\mathrm{AVE}}(x_{1},x_{2})=\eta_{\Delta}(x_{1},x_{2})\right]=0$.

Because $\phi_{\alpha}$ is independent of $(x_{1},x_{2})$ (note that the
distribution of $\epsilon^{\mathrm{AVE}}$ may depend on $(x_{1},x_{2})$ when
$\beta_{m}\neq 0$), we discuss the impact of $(x_{1},x_{2})$ on the OEP.
First, the OEP monotonically increases as the difference between the ground-
truth objective values, $\lvert f(x_{1};\Delta)-f(x_{2};\Delta)\rvert$,
increases. Second, the OEP is higher if
$\|\Gamma\nabla_{\bm{\epsilon}}f(x)\|_{\alpha}$ is lower for both $x_{1}$ and
$x_{2}$, where $\|\Gamma\nabla_{\bm{\epsilon}}f(x)\|_{\alpha}$ determines the
scaling factor for the randomness of the objective function. Indeed, if
$\alpha=2$, the scaling factor is the standard deviation of the objective
function value, $f(x;{\bm{\epsilon}})$. The OEP is higher if the ground-truth
objective function difference (i.e., signal) relative to noise scaling is
greater.

##  III Sign Averaging

A negative explicit averaging result occurs when the noise distribution has
stability parameter $\alpha\in(0,1]$. This is because expected value
$\mathbb{E}_{{\bm{\epsilon}}}[f(x;{\bm{\epsilon}})]$ does not exist, and the
average value does not converge. Hence, we propose an alternative strategy for
overcoming these limitations under different assumptions. The proposed _sign
averaging_ takes the average sign of the differences in the objective function
values of two candidate solutions. We theoretically reveal that sign averaging
allows to estimate the order of the medians of $f(x_{1};{\bm{\epsilon}})$ and
$f(x_{2};{\bm{\epsilon}})$ even when expected values
$\mathbb{E}[f(x_{1};{\bm{\epsilon}})]$ and
$\mathbb{E}[f(x_{2};{\bm{\epsilon}})]$ do not exist. We also propose an
approach for incorporating sign averaging into a comparison-based multipoint
search algorithm.

###  III-A Problem Statement

We describe the optimization objective and assumption on $f$ as follows.
Instead of assuming the linearity of $f$ with respect to ${\bm{\epsilon}}$ (1)
and a stable distribution of ${\bm{\epsilon}}$ (2), we make the following
technical assumption. The median difference between noisy objective function
values $f(x;{\bm{\epsilon}}_{x})$ and $f(y;{\bm{\epsilon}}_{y})$ is assumed to
be the difference between their medians.

######  Assumption 3 (Additivity of median).

Let ${\bm{\epsilon}}\in\mathbb{R}^{M}$ be a random vector and
$f(x;{\bm{\epsilon}}):\mathbb{R}^{D}\times\mathbb{R}^{M}\to\mathbb{R}$. In
addition, let $\phi_{f}(t;x)=\Pr_{\bm{\epsilon}}[f(x;{\bm{\epsilon}})\leqslant
t]$ be the cumulative distribution function of $f(x;{\bm{\epsilon}})$. Let
$m_{f}(x)$ be the center of the median of $f(x;{\bm{\epsilon}})$ defined as

$m_{f}(x)=\frac{\inf_{t}\\{\phi_{f}(t;x)>1/2\\}+\sup_{t}\\{\phi_{f}(t;x)<1/2\\}}{2}.$
(20)

Let $m_{f}(x,y)$ be the center of the median of
$f(x;{\bm{\epsilon}})-f(y;{\bm{\epsilon}}^{\prime})$, where
${\bm{\epsilon}}^{\prime}$ is an i.i.d. copy of $\epsilon$. Then,
$m_{f}(x,y)=m_{f}(x)-m_{f}(y)$.

In general, the median does not satisfy additivity
$m_{f}(x,y)=m_{f}(x)-m_{f}(y)$. However, 3 holds, for example, when the
distribution of $f(x;{\bm{\epsilon}})$ is symmetric for all
$x\in\mathbb{R}^{D}$, as formally stated below. The proof is provided in
Section A-D.

######  Proposition 5 (Sufficient condition for 3).

Suppose that the distribution of $f(x;{\bm{\epsilon}})$ is symmetric around
$m_{f}(x)$ for any $x\in\mathbb{R}^{D}$. Then, 3 holds.

3 holds also when the distribution of $f(x;{\bm{\epsilon}})-m_{f}(x)$ is
identical for all $x\in\mathbb{R}^{D}$. This case includes an often assumed
situation, namely the additive noise situation, where the noisy objective
function $f(x;{\bm{\epsilon}})=f(x)+\xi(\bm{\epsilon})$ is the sum of the
ground truth function $f(x)$ and the noise $\xi(\bm{\epsilon})$ that is
independent of $x$. It is formally stated below. The proof is provided in
Section A-E.

######  Proposition 6 (Sufficient condition for 3).

Suppose that the distribution of $f(x;{\bm{\epsilon}})-m_{f}(x)$ is identical
for all $x\in\mathbb{R}^{D}$. Then, 3 holds.

The median is not necessarily unique, but its center is, thus hindering
analysis. To mitigate this technical difficulty, we assume the following
properties.

######  Assumption 4 (Regularity of $\phi_{f}$).

For any $x\in\mathbb{R}^{D}$, the cumulative distribution function,
$\phi_{f}(t;x)$, is differentiable at $t=m_{f}(x)$, and its derivative (i.e.,
probability density function) is strictly positive at $t=m_{f}(x)$.

4 is not very restrictive. If 4 is not satisfied directly by $f(x;\epsilon)$,
it is satisfied by injecting an arbitrarily small noise to each query, as
formally stated below. The proof is provided in Section A-F.

######  Proposition 7 (Sufficient condition for 4).

If we replace each query with $x\mapsto f(x;\bm{\epsilon})+\xi$, where $\xi$
is an absolutely continuous random variable whose probability density is
continuous and strictly positive over $\mathbb{R}$, and $\xi$ and
$\bm{\epsilon}$ are independent, then 4 holds.

These two assumptions are satisfied for various symmetric distributions. For
example, under 1 and 2, the two assumptions above are satisfied when
$\beta_{m}=0$ for all $m=1,\dots,M$ (see Theorem 1.9 and Lemma 1.10 of [27]).
Moreover, these assumptions hold when neither 1 nor 2 is satisfied. Some
representative probability distributions satisfying 3 and 4 are the continuous
uniform, Student’s t, Laplace, Bates, and arcsine distributions as well as
part of the beta distribution.

In the following analysis, optimization index
$\mathcal{A}[f(x;{\bm{\epsilon}})]$ for minimization is the center of median
$m_{f}(x)$ of the objective function. In contrast to the optimization index
considered in Section II, $f(x;\Delta)$, which requires 2, and the expected
value, $\mathbb{E}_{\bm{\epsilon}}[f(x;{\bm{\epsilon}})]$, which requires the
existence of the expected value for all $x\in\mathbb{R}^{D}$, the median
always exists, and its center is uniquely determined. Although these concepts
differ, they are related. Under 1 and 2 with $\beta_{m}=0$ for all
$m=1,\dots,M$, we have
$f(x;\Delta)=m_{f}(x)=\mathbb{E}_{{\bm{\epsilon}}}[f(x;{\bm{\epsilon}})]$ for
all $\alpha\in(1,2]$ and $f(x;\Delta)=m_{f}(x)$ for all $\alpha\in(0,1]$,
whereas $\mathbb{E}_{\bm{\epsilon}}[f(x;{\bm{\epsilon}})]$ does not exist.

###  III-B OEP

Like in Section II, we investigate the OEP on the median as follows.

######  Definition 3 (OEP on median).

Let $\widehat{\eta}_{{\bm{\epsilon}}}(x_{1},x_{2})\in\\{-1,0,1\\}$ be an
estimate of

$\eta_{\mathrm{med}}(x_{1},x_{2}):=\mathrm{sign}\left(m_{f}(x_{1})-m_{f}(x_{2})\right).$
(21)

The _OEP_ on the median is given by

$\Pr_{\bm{\epsilon}}[\widehat{\eta}_{\bm{\epsilon}}(x_{1},x_{2})=\eta_{\mathrm{med}}(x_{1},x_{2})].$
(22)

We define the order estimate based on sign averaging as

$\widehat{\eta}_{\bm{\epsilon}}^{\mathrm{SA}}(x_{1},x_{2}):=\mathrm{sign}\left(\frac{1}{K}\sum_{k}^{K}\mathrm{sign}\left(f(x_{1};{\bm{\epsilon}}_{1,k})-f(x_{2};{\bm{\epsilon}}_{2,k})\right)\right).$
(23)

Consider its expected value:

$\mathbb{E}\left[\mathrm{sign}\left(f(x_{1};{\bm{\epsilon}}_{1})-f(x_{2};{\bm{\epsilon}}_{2})\right)\right].$
(24)

It exists as a finite value, and its sample average must converge to the
expected value as $K$ increases. Moreover, the following proposition
guarantees that the sign of expected value
$\mathbb{E}\left[\mathrm{sign}\left(f(x_{1};{\bm{\epsilon}}_{1})-f(x_{2};{\bm{\epsilon}}_{2})\right)\right]$
always agrees with $\eta_{\mathrm{med}}(x_{1},x_{2})$. The proof is provided
in Section A-G.

######  Proposition 8 (Expectation of
$\mathrm{sign}\left(f(x_{1};{\bm{\epsilon}}_{1})-f(x_{2};{\bm{\epsilon}}_{2})\right)$).

Let random vector ${\bm{\epsilon}}\in\mathbb{R}^{M}$ and objective function
$f(x;{\bm{\epsilon}})$ satisfy 3. Define

$\epsilon^{f}(x_{1},x_{2}):=f(x_{1};{\bm{\epsilon}}_{1})-f(x_{2};{\bm{\epsilon}}_{2}),$
(25)

where ${\bm{\epsilon}}_{1},{\bm{\epsilon}}_{2}$ are i.i.d. copies of
${\bm{\epsilon}}$. For any $(x_{1},x_{2})$, we have

$\mathbb{E}\left[\mathrm{sign}(\epsilon^{f}(x_{1},x_{2}))\right]\cdot\eta_{\mathrm{med}}(x_{1},x_{2})\geqslant
0.$ (26)

Moreover, if 4 is satisfied, we have the following relation for any
$(x_{1},x_{2})$ satisfying $\eta_{\mathrm{med}}(x_{1},x_{2})\neq 0$:

$\mathbb{E}\left[\mathrm{sign}(\epsilon^{f}(x_{1},x_{2}))\right]\cdot\eta_{\mathrm{med}}(x_{1},x_{2})>0.$
(27)

The following OEP bound with
$\widehat{\eta}_{\bm{\epsilon}}^{\mathrm{SA}}(x_{1},x_{2})$ guarantees that
the probability of
$\mathbb{E}\left[\mathrm{sign}(\epsilon^{f}(x_{1},x_{2}))\right]\cdot\eta_{\mathrm{med}}(x_{1},x_{2})>0$
being satisfied approaches $1$ as $K$ increases if
$\eta_{\mathrm{med}}(x_{1},x_{2})\neq 0$. The proof is presented in Section
A-H.

######  Theorem 9 (OEP with sign averaging).

Let random vector ${\bm{\epsilon}}\in\mathbb{R}^{M}$ and objective function
$f(x;{\bm{\epsilon}})$ satisfy 3. In addition, let
${\bm{\epsilon}}_{i,k}(i=1,2$ and $k=1,\ldots,K)$ be an i.i.d. copy of
${\bm{\epsilon}}$. If $\eta_{\mathrm{med}}(x_{1},x_{2})\neq 0$, it holds that

$\Pr_{\bm{\epsilon}}\left[\widehat{\eta}_{\bm{\epsilon}}^{\mathrm{SA}}(x_{1},x_{2})=\eta_{\mathrm{med}}(x_{1},x_{2})\right]\\\
\geqslant
1-\exp\left(-\frac{K\cdot\left(\mathbb{E}\left[\mathrm{sign}\left(\epsilon^{f}(x_{1},x_{2})\right)\right]\right)^{2}}{2}\right).$
(28)

The implications of Theorem 9 are as follows. First, we observe that the OEP
lower bound increases with increasing $K$. Because 3 and 4 are satisfied by 1
and 2 with $\beta_{m}=0$ (for all $m=1,\dots,M$), the result indicates that
sign averaging is effective even when $\alpha\leqslant 1$, whereas explicit
averaging fails. Hence, sign averaging is a valid alternative to explicit
averaging in comparison-based algorithms. Second, the OEP lower bound with
respect to pair $(x_{1},x_{2})$ is higher for a pair with larger
$\mathbb{E}\left[\mathrm{sign}\left(\epsilon^{f}(x_{1},x_{2})\right)\right]^{2}$.
This expectation can be rewritten as

$\mathbb{E}\left[\mathrm{sign}\left(\epsilon^{f}(x_{1},x_{2})\right)\right]\\\
=\Pr[f(x_{1};\epsilon_{1})<f(x_{2};\epsilon_{2})]-\Pr[f(x_{1};\epsilon_{1})>f(x_{2};\epsilon_{2})].$
(29)

Hence, a higher probability of a solution being better than or smaller than
the other one increases the OEP lower bound.

Considering Proposition 8, 4 guarantees (27) and the OEP lower bound with
$\eta_{\bm{\epsilon}}^{\mathrm{SA}}(x_{1},x_{2})$ can arbitrarily approach 1
by increasing $K$. Sample size $K$ for sign averaging to obtain an OEP above
$p$ is $O\left(\log\left(\frac{1}{1-p}\right)\right)$, as formally stated
below.

######  Corollary 10 (Sufficient sample size $K$ to estimate OEP).

Suppose that 3 and 4 hold. Then, for any pair $(x_{1},x_{2})$ satisfying
$\eta_{\mathrm{med}}(x_{1},x_{2})\neq 0$ and for any $p\in(0,1)$, we have
$\Pr_{\bm{\epsilon}}\left[\widehat{\eta}_{\bm{\epsilon}}^{\mathrm{SA}}(x_{1},x_{2})=\eta_{\mathrm{med}}(x_{1},x_{2}),\right]\geqslant
p$ if

$K\geqslant\frac{2}{\mathbb{E}\left[\mathrm{sign}\left(\epsilon^{f}(x_{1},x_{2})\right)\right]^{2}}\cdot\log\left(\frac{1}{1-p}\right).$
(30)

###  III-C Weighting Based on Sign Averaging

We propose a mechanism for incorporating sign averaging into comparison-based
optimization algorithms.

A comparison-based optimization algorithm typically has a set of candidate
solutions, $\\{x_{i}\\}_{i=1}^{\lambda}$, which is often called the
population. If objective function $f$ is deterministic, the population is
sorted in ascending order of $f$, that is,
$f(x_{1:\lambda})\leqslant\dots\leqslant f(x_{\lambda:\lambda})$, where
$i:\lambda$ denotes the index of the $i$-th best point among $\lambda$ points.
Ties (i.e., $f(x_{i})=f(x_{j})$) are usually ignored because they are unlikely
during continuous optimization. Based on the rankings, predefined weights are
assigned to the candidate solutions. Let $w_{1},\dots,w_{\lambda}$ be
predefined weights. Typically, $w_{1}\geqslant\dots\geqslant w_{\lambda}$.
Then, $x_{i:\lambda}$ receives $w_{i}$. The internal information of the
optimization algorithm is updated using
$\\{(x_{i:\lambda},w_{i})\\}_{i=1}^{\lambda}$. For example, in CMA-ES [23, 21,
20], the internal information comprises the parameters of the normal
distribution and evolution paths. The predefined weights can differ for
updates of different internal parameters, like for the latest CMA-ES [2]
(i.e., when active covariance matrix adaptation [1] is used).

One technical difficulty in incorporating sign averaging into comparison-based
optimization algorithms is that $\widehat{\eta}_{\bm{\epsilon}}^{\mathrm{SA}}$
is not necessarily transitive. That is, for three points $x_{1},x_{2},x_{3}$,
$\widehat{\eta}_{\bm{\epsilon}}^{\mathrm{SA}}(x_{1},x_{2})<0$, and
$\widehat{\eta}_{\bm{\epsilon}}^{\mathrm{SA}}(x_{2},x_{3})<0$,
$\widehat{\eta}_{\bm{\epsilon}}^{\mathrm{SA}}(x_{1},x_{3})>0$ may occur. In
such case, the rankings of $x_{1},x_{2},x_{3}$ cannot be determined. If this
occurs, weighting cannot be directly applied. To sort population
$\\{x_{i}\\}_{i=1}^{\lambda}$ using
$\widehat{\eta}_{\bm{\epsilon}}^{\mathrm{SA}}$, the score of each point
$x_{i}$ is defined as the number of points that are estimated to be better
than $x_{i}$, namely,

$R\left(x_{i};\\{x_{j}\\}_{j=1}^{\lambda}\right):=\sum_{j=1}^{\lambda}\mathbbm{1}{\left\\{\widehat{\eta}_{\bm{\epsilon}}^{\mathrm{SA}}(x_{j},x_{i})\leqslant
0\right\\}},$ (31)

and the population can be sorted in ascending order of $R$. If the population
rankings are uniquely determined by
$\widehat{\eta}_{\bm{\epsilon}}^{\mathrm{SA}}$ (i.e.,
$\widehat{\eta}_{\bm{\epsilon}}^{\mathrm{SA}}$ is transitive on
$\\{x_{j}\\}_{j=1}^{\lambda}$), (31) is the ranking of $x_{i}$. Therefore, it
may be considered as an extension of the ranking mechanism using
$\widehat{\eta}_{\bm{\epsilon}}^{\mathrm{SA}}$. However, ties may occur and
should be handled carefully to avoid undermining the performance of the
resulting algorithm.

In [5], a weighting scheme was proposed to treat ties with mathematical rigor
for deterministic objective functions. Weighting introduced in [5] is defined
as follows. For each $x_{i}\in\\{x_{j}\\}_{j=1}^{\lambda}$, the numbers of
strictly better and weakly better points among $\\{x_{j}\\}_{j=1}^{\lambda}$
are counted as follows:

$\displaystyle
r^{<}_{i}:=\sum_{j=1}^{\lambda}\mathbbm{1}{\left\\{f(x_{j})<f(x_{i})\right\\}},$
(32) $\displaystyle
r^{\leqslant}_{i}:=\sum_{j=1}^{\lambda}\mathbbm{1}{\left\\{f(x_{i})\leqslant
f(x_{j})\right\\}}.$ (33)

Based on these quantities, a weight is assigned to every $x_{i}$ as follows:

$\bar{w}_{i}:=\frac{1}{r^{\leqslant}_{i}-r^{<}_{i}}\cdot\sum_{j=r^{<}_{i}+1}^{r^{\leqslant}_{i}}w_{j}.$
(34)

If there is no tie, $r_{i}^{\leqslant}=r_{i}^{<}+1$ for all
$i=1,\dots,\lambda$ and $r_{i}^{\leqslant}$ agrees with the conventional
$x_{i}$ ranking. The assigned weights are equivalent to the aforementioned
conventional weighting. If ties exist, the same weight is assigned to the tie
points. Thus, we can guarantee that the sum of assigned weights is the same as
the sum of the predefined weights, that is,
$\sum_{i=1}^{\lambda}\bar{w}_{i}=\sum_{i=1}^{\lambda}w_{i}$.

We incorporate averaging into weighting in (34) by simply replacing $f(\cdot)$
with $R\left(\cdot;\\{x_{j}\\}_{j=1}^{\lambda}\right)$ in (32) and (33). If
$f$ is deterministic and no ties exist among $f(x_{1}),\dots,f(x_{\lambda})$,
the assigned weights are the same as those obtained by conventional weighting.
If $f$ is stochastic but the rankings of the population are uniquely
determined by $\widehat{\eta}_{\bm{\epsilon}}^{\mathrm{SA}}$ (i.e.,
$\widehat{\eta}_{\bm{\epsilon}}^{\mathrm{SA}}$ is transitive on
$\\{x_{j}\\}_{j=1}^{\lambda}$), the weights computed by the proposed approach
are the same as those assigned by conventional weighting with
$\widehat{\eta}_{\bm{\epsilon}}^{\mathrm{SA}}$ for comparison. Moreover,
according to Theorem 9, the weights assigned by the proposed approach agree
with conventional weighting with $\eta_{\mathrm{med}}$ for comparison in limit
$K\to\infty$ if no tie exists.

##  IV Experiments on OEP

Figure 1: Results of 10 runs of CMA-ES with explicit averaging. The solid
lines indicate the median of the moving average of Tau-b (left axis) with
sample sizes $K=1$ ($\bullet$), $K=10$ ($\blacktriangledown$), and $K=50$
($\blacksquare$). The span of the moving average is $10$. The dashed lines
indicate $-\log(f(m_{t};\Delta))$ (right axis) with sample sizes $K=1$
($\star$), $K=10$ ($\times$), and $K=50$ ($\blacklozenge$). The top and bottom
of the band correspond to the 75% and 25% values among trials. The medians are
taken over $10$ trials. (ANE, additive-noise ellipsoid; LNE, linear-noise
ellipsoid; MNE, multiplicative-noise ellipsoid)  Figure 2: Results of 10 runs
of CMA-ES with sign averaging. The solid lines indicate the median of the
moving average of Tau-b (left axis) with sample sizes $K=1$ ($\bullet$),
$K=10$ ($\blacktriangledown$), and $K=50$ ($\blacksquare$). The dashed lines
indicate $-\log(f(m_{t};\Delta))$ (right axis) with sample sizes $K=1$
($\star$), $K=10$ ($\times$), and $K=50$ ($\blacklozenge$).

First, we incorporated explicit averaging into $(\mu/\mu_{w},\lambda)$-CMA-ES
[23, 21] as a representative comparison-based optimization algorithm. The
negative effect of averaging was confirmed when the value of stability
parameter $\alpha$ was less than $1$ and positive for $\alpha\in(1,2]$. Under
the same experimental settings, explicit averaging was replaced by sign
averaging using the weighting scheme discussed in Section III-C. Sign
averaging contributed to the CMA-ES accuracy for approximating the optimum,
even for $\alpha\in(0,1]$ at which conventional averaging failed.

In the abovementioned experiments, two optimization indices were obtained per
iteration. The first index was the objective function value with location
parameter $f(m_{t},\Delta)$, which was assumed to be the ground-truth $f$
value. Again, this index agreed with both the median of $f(m_{t},\epsilon)$ on
symmetric stable noise and the mean if it was finite. The second index was the
Kendall rank correlation coefficient, Tau-b, which is defined on pairs of
numerical values as the difference between the number of concordant pairs and
the number of discordant pairs divided by the number of pairs that are not
tied [3]. In each experiment, Tau-b was applied to a set of pairs of the
estimated ranking of the ground-truth $f$ value, $f(x;\Delta)$, and the proper
ranking in the population at each iteration was calculated. Tau-b is 1 if
noise handling attains perfect ranking estimation and $-$1 otherwise, while it
is 0 if the estimated and ground-truth rankings are uncorrelated.

#### Test functions

We define three test functions that satisfy 1. 1) Additive-noise ellipsoid
$f_{\mathrm{add}}(x;{\bm{\epsilon}})$ is defined as
$f_{\mathrm{add}}(x;{\bm{\epsilon}})=x^{\mathrm{T}}Hx+\epsilon$, where
$H\in\mathbb{R}^{D\times D}$ is a positive definite diagonal matrix with the
$i$-th eigenvalue being $10^{2(i-1)/(D-1)}$, and $\epsilon\in\mathbb{R}$ is a
random variable following a stable distribution. The values of these
parameters are listed in the titles of the corresponding figures.

2) Multiplicative-noise ellipsoid $f_{\mathrm{mul}}(x;{\bm{\epsilon}})$ is
defined as
$f_{\mathrm{mul}}(x;{\bm{\epsilon}})=x^{\mathrm{T}}Hx\cdot\epsilon$, where
$H\in\mathbb{R}^{D\times D}$ is identical to that of the additive-noise
ellipsoid and $\epsilon\in\mathbb{R}$ is stable noise.

3) Linear-noise ellipsoid $f_{\mathrm{LNE}}(x;{\bm{\epsilon}})$ is defined as
$f_{\mathrm{LNE}}(x;{\bm{\epsilon}})=x^{\mathrm{T}}Hx+{\bm{\epsilon}}^{\mathrm{T}}x$,
where $H$ is identical to that of the additive noise ellipsoid and
${\bm{\epsilon}}\in\mathbb{R}^{M}$ is a random vector whose $m$-th element is
stable noise ${\bf S}(\alpha,\beta,\gamma_{m},\delta_{m})$. We set
$\gamma_{m}=10^{a+(m-1)(b-a)/(D-1)}$ and $\delta_{m}=0$. The values of $a,b$
are specified as ${\bm{S}}(\alpha,\beta,[10^{a},10^{b}],[0,0])$ in the titles
of the corresponding figures.

For the three defined test functions, the ground-truth objective function
value is $f(x;\Delta)=x^{\mathrm{T}}Hx$ as we set $\delta_{m}=0$ for all
$m=1,\dots,M$.

#### Parameters of noise ${\bm{\epsilon}}$

We investigated the ability of algorithms to optimize the three test functions
under five types of stable noise components: 1) normal distribution
($\alpha=2,\beta$), 2) symmetric stable distribution with
$\alpha=1.5,\beta=0$, 3) Cauchy distribution ($\alpha=1,\beta=0$), 4)
symmetric stable distribution with $\alpha=0.5,\beta=0$, and 5) Lévy
distribution ($\alpha=0.5,\beta=1$). The first four distributions were
examined to compare the differences between variations in $\alpha$. The Lévy
distribution was used to evaluate an asymmetric stable distribution.

#### Other experimental settings

For each test function, the dimension of design variable $x$ was $D=20$.
Initial solution $m_{0}$ was $[10,\ldots,10]^{\mathrm{T}}$ for all the runs.
All the settings of CMA–ES followed the details in [20].

###  IV-A Explicit Averaging

Averaging was incorporated into CMA-ES. It received
$\\{\frac{1}{K}\sum_{k}^{K}f(x_{i};\epsilon_{i,k})\\}_{i}^{\lambda}$ as the
estimated $f$ value of population $\\{x_{i}\\}_{i}^{\lambda}$. Figure 1 shows
Tau-b on pairs
$\left\\{\left(\frac{1}{K}f(x_{i};\epsilon_{i,k}),f(x_{i};\Delta)\right)\right\\}_{i}^{\lambda}$
and ground-truth objective function value $-\log(f(m_{t};\Delta)$ over
iteration $t$ of the CMA-ES with explicit averaging for $K=1,10,50$, where
$m_{t}$ is the mean vector of CMA-ES at iteration $t$. All the test functions
and noise satisfied 1 and 2.

In Figure 1, the implication of Theorem 3 on sample size $K$ was confirmed. 1)
For $\alpha\in(1,2]$, we obtained a lower $f(x;\Delta)$ for a greater $K$
(first and second columns from the left). 2) For $\alpha=1$, averaging had no
impact on $\tau$ and $f(x;\Delta)$, and it only increased the number of $f$
calls (third column). 3) For $\alpha\in(0,1)$, we obtained a lower
$f(x;\Delta)$ for a lower $K$ (fourth and fifth columns).

###  IV-B Sign Averaging

Next, we incorporated sign averaging into CMA-ES using the weighting proposed
in Section III-C, obtaining sign-averaging CMA-ES. We confirmed the
effectiveness of sign averaging for problems with noise without a finite mean,
for which explicit averaging was unsuccessful. For a fair comparison, all the
experimental settings and baseline algorithm were as those in Section IV,
except for explicit averaging, which was replaced by sign averaging. Because
sign averaging did not directly estimate the $f$ values, Tau-b on the pairs of
the estimated rankings of $x_{i}$ in the population, $R(x_{i})$, and the
ground-truth value were examined. All the test functions with symmetric
distributions of the noise used in Section IV also satisfied 3 and 4, and
asymmetric noise ($\alpha=0.5,\beta=1$) did not satisfy 3.

Figure 2 shows Kendall’s Tau-b on
$\left(R(x_{i}),f(x_{i};\Delta)\right)_{i}^{\lambda}$ and
$-\log(f(m_{t};\Delta))$ over iteration $t$. For $\alpha=2$ (first column in
Figure 2), sign-averaging CMA-ES attained almost the same precision of
$-\log(f(m_{t};\Delta))$ as that obtained by CMA-ES with explicit averaging
(Figure 1) for all $K$ values in the three test functions. For $\alpha=1.5$
(second column in Figure 2), the performance of sign-averaging CMA-ES remained
almost the same as that for $\alpha=2$ regardless of the $K$ value and test
function. For the other two symmetric distributions (third and fourth columns
in Figure 2), sign-averaging CMA-ES clearly outperformed CMA-ES with explicit
averaging on the three test functions as $K$ increased. This finding agreed
with the implications of Theorem 9. Furthermore, for asymmetric noise (last
column in Figure 2), while no effect of $K$ appeared in the multiplicative-
noise ellipsoid, the performance on index $f(x;\Delta)$ was improved by
increasing $K$ in the additive- and linear-noise ellipsoids.

##  V Discussion

###  V-A Invariance under Strictly Increasing Transformation

Let $G:\mathbb{R}\to\mathbb{R}$ be a strictly increasing transformation such
that $G(x)\lesseqgtr G(y)\Leftrightarrow x\lesseqgtr y$. The behavior of a
comparison-based algorithm is generally invariant to such transformation
because the ranking of the objective function values is invariant to $G$.
However, invariance may break with averaging, which is not invariant to $G$.
Hence, the following case may occur:

$\mathrm{sign}\left\\{\frac{1}{K}\sum_{k}^{K}G\circ
f(x_{1};{\bm{\epsilon}}_{1})-\frac{1}{K}\sum_{k}^{K}G\circ
f(x_{2};{\bm{\epsilon}}_{2})\right\\}\\\
\neq\mathrm{sign}\left\\{\frac{1}{K}\sum_{k}^{K}f(x_{1};{\bm{\epsilon}}_{1})-\frac{1}{K}\sum_{k}^{K}f(x_{2};{\bm{\epsilon}}_{2})\right\\}.$
(35)

On the other hand, sign averaging preserves invariance under $G$. Hence, the
following equality always holds:

$\mathrm{sign}\left\\{\frac{1}{K}\sum_{k}^{K}\mathrm{sign}\left(G\circ
f(x_{1};{\bm{\epsilon}}_{1})-G\circ
f(x_{2};{\bm{\epsilon}}_{2})\right)\right\\}\\\
=\mathrm{sign}\left\\{\frac{1}{K}\sum_{k}^{K}\mathrm{sign}\left(f(x_{1};{\bm{\epsilon}}_{1})-f(x_{2};{\bm{\epsilon}}_{2})\right)\right\\}.$
(36)

###  V-B Asymptotic OEP

Considering asymptotic normality, we can derive the order of the OEP in $K$
for explicit averaging, sign averaging, and sample median.

For a global perspective of the derivation, suppose that
$\widehat{\theta}_{K}$ is an estimate of a parameter $\theta$ with $K$ samples
and that $\widehat{\theta}_{K}$ exhibits asymptotic normality, that is,
$Z_{K}=\frac{\sqrt{K}(\widehat{\theta}_{K}-\theta)}{\sigma}$ converges in
distribution to $\mathcal{N}(0,1)$ for some $\sigma>0$. Then, for $\theta\neq
0$, we have

$\displaystyle\Pr\left[\mathrm{sign}(\widehat{\theta}_{K})=\mathrm{sign}(\theta)\right]$
$\displaystyle=\Pr\left[\widehat{\theta}_{K}\cdot\theta>0\right]$ (37a)
$\displaystyle=\Pr\left[\left(\theta+\frac{\sigma
Z_{K}}{\sqrt{K}}\right)\cdot\theta>0\right]$ (37b)
$\displaystyle=\Pr\left[Z_{K}\leqslant\frac{\sqrt{K}\lvert\theta\rvert}{\sigma}\right]$
(37c)
$\displaystyle\approx\Phi_{\mathcal{N}}\left(\frac{\sqrt{K}\lvert\theta\rvert}{\sigma}\right).$
(37d)

Therefore, the probability of correctly estimating the sign of parameter
$\theta$ is approximated by the cumulative distribution function of a normal
distribution.

The OEPs for explicit averaging, sign averaging, and sample median are
approximated using asymptotic normality. For simplicity, we assume that the
distribution of $f(x;{\bm{\epsilon}})$ is symmetric and absolutely continuous
for any $x$ and that the probability density function of
$f(x_{1};{\bm{\epsilon}}_{1})-f(x_{2};{\bm{\epsilon}}_{2})$ is strictly
positive at its unique median. Hence, the following results are obtained:

  * •

(Explicit averaging) Let
$\widehat{\mu}_{K}^{\mathrm{AVE}}(x_{1},x_{2}):=\frac{1}{K}\sum_{k}^{K}f(x_{1};{\bm{\epsilon}}_{1,k})-\frac{1}{K}\sum_{k}^{K}f(x_{2};{\bm{\epsilon}}_{2,k})$,
$\mu^{\mathrm{AVE}}(x_{1},x_{2}):=\mathbb{E}[f(x_{1};{\bm{\epsilon}}_{1})-f(x_{2};{\bm{\epsilon}}_{2})]$,
and
$\sigma^{\mathrm{AVE}}(x_{1},x_{2})^{2}:=\operatorname{Var}[f(x_{1};{\bm{\epsilon}})]+\operatorname{Var}[f(x_{2};{\bm{\epsilon}})]$.
Assume $\mu^{\mathrm{AVE}}(x_{1},x_{2})<\infty$ and
$\sigma^{\mathrm{AVE}}(x_{1},x_{2})<\infty$ for any $x$. Then, considering the
central limit theorem, the asymptotic normality of
$\widehat{\mu}^{\mathrm{AVE}}(x_{1},x_{2})$ for
$\mu^{\mathrm{AVE}}(x_{1},x_{2})$ and $\sigma^{\mathrm{AVE}}(x_{1},x_{2})$
holds, and for large $K$, the OEP is approximated as

$\Pr\left[\mathrm{sign}(\widehat{\mu}^{\mathrm{AVE}}(x_{1},x_{2}))=\mathrm{sign}(\mu^{\mathrm{AVE}}(x_{1},x_{2}))\right]\\\
\approx\phi_{\mathcal{N}}\left(\frac{\sqrt{K}\cdot\lvert\mu^{\mathrm{AVE}}(x_{1},x_{2})\rvert}{\sigma^{\mathrm{AVE}}(x_{1},x_{2})}\right).$
(38)

  * •

(Sign averaging) Let
$\widehat{\mu}_{K}^{\mathrm{SA}}(x_{1},x_{2}):=\frac{1}{K}\sum_{k}^{K}\mathrm{sign}(f(x_{1};{\bm{\epsilon}}_{1,k})-f(x_{2};{\bm{\epsilon}}_{2,k}))$,
$\mu^{\mathrm{SA}}:=\mathbb{E}[\mathrm{sign}(f(x_{1};{\bm{\epsilon}}_{1,k})-f(x_{2};{\bm{\epsilon}}_{2,k}))]$,
and
$\sigma^{\mathrm{SA}}(x_{1},x_{2})^{2}:=\operatorname{Var}[\mathrm{sign}(f(x_{1};{\bm{\epsilon}}_{1,k})-f(x_{2};{\bm{\epsilon}}_{2,k}))]$.
Then, $\mu^{\mathrm{SA}}(x_{1},x_{2})<\infty$ and
$\sigma^{\mathrm{SA}}(x_{1},x_{2})^{2}<\infty$. Considering the central limit
theorem, the asymptotic normality of
$\widehat{\mu}^{\mathrm{SA}}(x_{1},x_{2})$ for
$\mu^{\mathrm{SA}}(x_{1},x_{2})$ and $\sigma^{\mathrm{SA}}(x_{1},x_{2})$
holds. For large $K$, the OEP is approximated as

$\Pr\left[\mathrm{sign}(\widehat{\mu}^{\mathrm{SA}}(x_{1},x_{2}))=\mathrm{sign}(\mu^{\mathrm{SA}}(x_{1},x_{2}))\right]\\\
\approx\phi_{\mathcal{N}}\left(\frac{\sqrt{K}\cdot\lvert\mu^{\mathrm{SA}}(x_{1},x_{2})\rvert}{\sigma^{\mathrm{SA}}(x_{1},x_{2})}\right).$
(39)

  * •

(Sample median) Let $\xi(x_{1},x_{2})$ be the value of the probability density
function of $f(x_{1};{\bm{\epsilon}}_{1})-f(x_{2};{\bm{\epsilon}}_{2})$ at its
median. In addition, let $\widehat{\mu}^{\mathrm{MED}}_{K}(x_{1},x_{2})$ be
the sample median of
$\\{f(x_{1};{\bm{\epsilon}}_{1,k})-f(x_{2};{\bm{\epsilon}}_{2,k})\\}_{k=1}^{K}$,
$\mu^{\mathrm{MED}}(x_{1},x_{2})$ be the population median of
$f(x_{1};{\bm{\epsilon}}_{1})-f(x_{2};{\bm{\epsilon}}_{2})$, and
$\sigma^{\mathrm{MED}}(x_{1},x_{2}):=\frac{1}{2\xi(x_{1},x_{2})}$. Consider
the asymptotic normality of the sample median [13]. For large $K$, the OEP is
approximated as

$\Pr\left[\mathrm{sign}(\widehat{\mu}_{K}^{\mathrm{MED}}(x_{1},x_{2}))=\mathrm{sign}(\mu^{\mathrm{MED}}(x_{1},x_{2}))\right]\\\
\approx\phi_{\mathcal{N}}\left(\frac{\sqrt{K}\cdot\lvert\mu^{\mathrm{MED}}(x_{1},x_{2})\rvert}{\sigma^{\mathrm{MED}}(x_{1},x_{2})}\right).$
(40)

Note that $\mathrm{sign}(\mu^{\mathrm{SA}}(x_{1},x_{2}))$ and
$\mathrm{sign}(\mu^{\mathrm{MED}}(x_{1},x_{2}))$ are equivalent for any
$x_{1},x_{2}$ considering Proposition 8. In addition, they are equivalent to
$\mathrm{sign}(\mu^{\mathrm{AVE}}(x_{1},x_{2}))$ if
$\mu^{\mathrm{AVE}}(x_{1},x_{2})<\infty$. Therefore, the three methods
estimate the same quantity under the current assumptions. Overall, we have the
same asymptotic order of OEP with respect to $K$ for the three estimates if
the mean of the noise exists as a finite value. We emphasize that the sample
median and the sign averaging attains these asymptotic results even if the
mean is not defined.

##  VI Conclusion

We investigated the effects of explicit averaging on comparison-based search
algorithms. In particular, we analyzed the OEP (Definition 2) of correctly
estimating the order of the ground-truth objective functions of pairs of
points. Assuming linearity of the objective function with respect to the noise
vector (1) and stably distributed noise (2), the main results (Theorem 3 and
Corollary 4) reveal the following. 1) If the stability parameter of stably
distributed noise is $\alpha\in(1,2]$, increasing sample size $K$ in explicit
averaging increases the probability. 2) If $\alpha=1$ (e.g., Cauchy
distribution), $K$ has no effect on the probability. 3) If $0<\alpha<1$ (e.g.,
Lévy distribution), increasing $K$ decreases the probability. Thus, explicit
averaging is effective even if the variance of $f(x;{\bm{\epsilon}})$ does not
exist, but it is harmful if the mean of $f(x;{\bm{\epsilon}})$ does not exist.
This is experimentally confirmed through optimization using CMA-ES.

We propose the alternative of sign averaging for estimating the order of the
objective function medians on a pair of points even when the objective
function mean does not exist. This is theoretically demonstrated in Theorem 9
and Corollary 10 under certain regularity conditions of the noise distribution
(3 and 4). Any problem satisfying 1 and 2 also satisfies 3 and 4 if the
distribution of the noise is symmetric. That is, the sign averaging can handle
the stable noise for all $\alpha\in(0,2]$, which includes the cases that the
mean does not exists as a finite value, where the averaging fails. Because
sign averaging cannot be directly incorporated into common comparison-based
algorithms, we propose a versatile weighting scheme. The efficiency of sign
averaging in this weighting scheme is experimentally confirmed using the CMA-
ES.

In future work, we will address the limitations of this study. First, we will
demonstrate how sign-averaging CMA-ES works under diverse noise distributions.
In this work, the analysis of sign averaging is based on 3. It is satisfied
when the distribution of the noisy objective function values is symmetric or
when the noise distribution is independent of design variable $x$. These
assumptions may not hold in some cases. We have not investigated the effect of
sign averaging in situations where 3 does not hold. We will investigate the
effect of sign averaging both theoretically and empirically. Second, we will
theoretically and empirically investigate the effect of the sample median as
an alternative noise-handling technique. The sample median may be the most
natural choice for noise handling when the optimization index is the
population median of the objective function. However, we have not investigated
the effect of the sample median under a finite sample size $K$. We will
theoretically and empirically investigate the effect of the sample median and
compare it with sign averaging.

##  Appendix A Proof

###  A-A Proof of Lemma 2

###### Proof.

First, we prove that the distribution of the sum of $\epsilon_{l}\sim{\bf
S}(\alpha,\beta_{l},\gamma_{l},\delta_{l})(l=1,\ldots,L)$ is given by

$\sum_{l}^{L}\epsilon_{l}\sim{\bf
S}\left(\alpha,\frac{\sum_{l}^{L}\beta_{l}\gamma_{l}^{\alpha}}{\sum_{l}^{L}\gamma_{l}^{\alpha}},\left(\sum_{l}^{L}\gamma_{l}^{\alpha}\right)^{\frac{1}{\alpha}},\sum_{l}^{L}\delta_{l}\right).$
(41)

For $L=2$, (41) is directly derived from Proposition 1. Assume that (41) holds
for $\sum_{l}^{L-1}\epsilon_{l}$. Then, considering Proposition 1, the scale
parameter of $\sum_{l}^{L-1}\epsilon_{l}+\epsilon_{L}$ is given by
$(\sum_{l}^{L-1}\gamma_{l}^{\alpha}+\gamma_{L}^{\alpha})^{\frac{1}{\alpha}}=(\sum_{l}^{L}\gamma_{l}^{\alpha})^{\frac{1}{\alpha}}$.
The skewness parameter is given by

$\displaystyle\frac{\frac{\sum_{l}^{L-1}\beta_{l}\gamma_{l}^{\alpha}}{\sum_{l}^{L-1}\gamma_{l}^{\alpha}}\cdot\sum_{l}^{L-1}\gamma_{l}^{\alpha}+\frac{\beta_{L}\gamma_{L}^{\alpha}}{\gamma_{L}^{\alpha}}\cdot\gamma_{L}^{\alpha}}{\sum_{l}^{L}\gamma_{l}^{\alpha}}$
$\displaystyle=\frac{\sum_{l}^{L-1}\beta_{l}\gamma_{l}^{\alpha}+\beta_{L}\gamma_{L}^{\alpha}}{\sum_{l}^{L}\gamma_{l}^{\alpha}}$
(42a)
$\displaystyle=\frac{\sum_{l}^{L}\beta_{l}\gamma_{l}^{\alpha}}{\sum_{l}^{L}\gamma_{l}^{\alpha}}.$
(42b)

The location parameter is given by
$\sum_{l}^{L-1}\delta_{l}+\delta_{L}=\sum_{l}^{L}\delta_{l}$. This completes
the proof of (41).

Considering Taylor’s theorem, because $f$ is linear over ${\bm{\epsilon}}$, we
have ’

$f(x_{1};{\bm{\epsilon}}_{1})-f(x_{2};{\bm{\epsilon}}_{2})=f(x_{1};\Delta)+\nabla_{\bm{\epsilon}}f(x_{1})^{\mathrm{T}}({\bm{\epsilon}}_{1,k}-\Delta)\\\
-f(x_{2};\Delta)-\nabla_{\bm{\epsilon}}f(x_{2})^{\mathrm{T}}({\bm{\epsilon}}_{2,k}-\Delta).$
(43)

Consider the distribution of
$\nabla_{\bm{\epsilon}}f(x)^{\mathrm{T}}({\bm{\epsilon}}_{1,k}-\Delta)$. From
Proposition 1,

$g_{m}(x)\cdot(\epsilon_{m}-\delta_{m})\sim{\bf
S}\Big{(}\alpha,\mathrm{sign}(g_{m}(x))\cdot\beta_{m},|g_{m}(x)|\cdot\gamma_{m},\\\
-\frac{2}{\pi}\beta_{m}\gamma_{m}g_{m}(x)\log|g_{m}(x)|\cdot\mathbbm{1}{\left\\{\alpha=1\right\\}}\Big{)}.$
(44)

Hence, using (41),

$\nabla_{\bm{\epsilon}}f(x)^{\mathrm{T}}({\bm{\epsilon}}-\Delta)=\sum_{m}^{M}g_{m}(x)\cdot(\epsilon_{m}-\delta_{m})\\\
\sim{\bf
S}\left(\alpha,\beta^{\prime}(x),\gamma^{\prime}(x),\delta(x)^{\prime}\cdot\mathbbm{1}{\left\\{\alpha=1\right\\}}\right)\enspace.$
(45)

Moreover, according to Proposition 1,

$\nabla_{\bm{\epsilon}}f(x_{1})^{\mathrm{T}}({\bm{\epsilon}}_{1}-\Delta)-\nabla_{\bm{\epsilon}}f(x_{2})^{\mathrm{T}}({\bm{\epsilon}}_{2}-\Delta)\\\
\sim{\bf
S}\left(\alpha,\beta^{\prime\prime}(x_{1},x_{2}),\gamma^{\prime\prime}(x_{1},x_{2}),\delta^{\prime\prime}(x_{1},x_{2})\cdot\mathbbm{1}{\left\\{\alpha=1\right\\}}\right)\enspace.$
(46)

Then, by adding $f(x_{1};\Delta)-f(x_{2};\Delta)$, we obtain

$f(x_{1};{\bm{\epsilon}}_{1})-f(x_{2};{\bm{\epsilon}}_{2})\sim{\bf
S}\Big{(}\alpha,\beta^{\prime\prime}(x_{1},x_{2}),\gamma^{\prime\prime}(x_{1},x_{2})\\\
,f(x_{1};\Delta)-f(x_{2};\Delta)+\delta^{\prime\prime}(x_{1},x_{2})\cdot\mathbbm{1}{\left\\{\alpha=1\right\\}}\Big{)}.$
(47)

Again, Proposition 1 and (41) suggest that

$\frac{1}{K}\sum_{k}^{K}\left(f(x_{1};{\bm{\epsilon}}_{1})-f(x_{2};{\bm{\epsilon}}_{2})\right)\\\
\sim{\bf
S}\Big{(}\alpha,\beta^{\prime\prime}(x_{1},x_{2}),K^{\frac{1}{\alpha}-1}\cdot\gamma^{\prime\prime}(x_{1},x_{2}),\\\
f(x_{1};\Delta)-f(x_{2};\Delta)+\delta^{\prime\prime}(x_{1},x_{2})\cdot\mathbbm{1}{\left\\{\alpha=1\right\\}}\Big{)}.$
(48)

This completes the proof. ∎

###  A-B Proof of Theorem 3

###### Proof.

From Lemma 2 and Proposition 1, we have

$\frac{1}{K^{\frac{1}{\alpha}-1}\cdot\gamma^{\prime\prime}(x_{1},x_{2})}\cdot\left(\frac{1}{K}\sum_{k}^{K}f(x_{1};{\bm{\epsilon}}_{1,k})-\frac{1}{K}\sum_{k}^{K}f(x_{2};{\bm{\epsilon}}_{2,k})\right)\\\
\sim{\bm{S}}\Bigg{(}\alpha,\beta^{\prime\prime}(x_{1},x_{2}),1,\frac{f(x_{1};\Delta)-f(x_{2};\Delta)}{K^{\frac{1}{\alpha}-1}\cdot\gamma^{\prime\prime}(x_{1},x_{2})}\\\
+\left(\frac{\delta^{\prime\prime}(x_{1},x_{2})}{\gamma^{\prime\prime}(x_{1},x_{2})}+\frac{2}{\pi}\beta^{\prime\prime}(x_{1},x_{2})\cdot\log(\gamma^{\prime\prime}(x_{1},x_{2}))\right)\cdot\mathbbm{1}{\left\\{\alpha=1\right\\}}\Bigg{)}.$
(49)

$K^{\frac{1}{\alpha}-1}$ disappears in the second item of the location
parameter for $\alpha=1$ and $\gamma^{\prime\prime}(x_{1},x_{2})>0$ for any
$x_{1},x_{2}$. As
$\epsilon^{\mathrm{AVE}}+\frac{f(x_{1};\Delta)-f(x_{2};\Delta)}{K^{\frac{1}{\alpha}-1}\cdot\gamma^{\prime\prime}(x_{1},x_{2})}$
and
$\frac{1}{K^{\frac{1}{\alpha}-1}\cdot\gamma^{\prime\prime}(x_{1},x_{2})}\cdot\left(\frac{1}{K}\sum_{k}^{K}f(x_{1};{\bm{\epsilon}}_{1,k})-\frac{1}{K}\sum_{k}^{K}f(x_{2};{\bm{\epsilon}}_{2,k})\right)$
are equivalent in distribution, we have

$\displaystyle\Pr_{\bm{\epsilon}}\left[\widehat{\eta}_{\epsilon}^{\mathrm{AVE}}(x_{1},x_{2})=\eta_{\Delta}(x_{1},x_{2})\right]$
(50a) $\displaystyle=\begin{aligned}
\Pr_{\bm{\epsilon}}\Bigg{[}&\left(\frac{1}{K}\sum_{k}^{K}\left(f(x_{1},{\bm{\epsilon}}_{1,k})-f(x_{2},{\bm{\epsilon}}_{2,k})\right)\right)\\\
&\cdot(f(x_{1};\Delta)-f(x_{2};\Delta))>0\Bigg{]}\end{aligned}$ (50b)
$\displaystyle\begin{multlined}=\Pr_{\bm{\epsilon}}\left[\epsilon^{\mathrm{AVE}}>-\frac{f(x_{1};\Delta)-f(x_{2};\Delta)}{K^{\frac{1}{\alpha}-1}\cdot\gamma^{\prime\prime}(x_{1},x_{2})}\right]\cdot\mathbbm{1}{\left\\{\eta_{\Delta}(x_{1},x_{2})>0\right\\}}\\\
+\Pr_{\bm{\epsilon}}\left[\epsilon^{\mathrm{AVE}}<-\frac{f(x_{1};\Delta)-f(x_{2};\Delta)}{K^{\frac{1}{\alpha}-1}\cdot\gamma^{\prime\prime}(x_{1},x_{2})}\right]\cdot\mathbbm{1}{\left\\{\eta_{\Delta}(x_{1},x_{2})<0\right\\}}\end{multlined}=\Pr_{\bm{\epsilon}}\left[\epsilon^{\mathrm{AVE}}>-\frac{f(x_{1};\Delta)-f(x_{2};\Delta)}{K^{\frac{1}{\alpha}-1}\cdot\gamma^{\prime\prime}(x_{1},x_{2})}\right]\cdot\mathbbm{1}{\left\\{\eta_{\Delta}(x_{1},x_{2})>0\right\\}}\\\
+\Pr_{\bm{\epsilon}}\left[\epsilon^{\mathrm{AVE}}<-\frac{f(x_{1};\Delta)-f(x_{2};\Delta)}{K^{\frac{1}{\alpha}-1}\cdot\gamma^{\prime\prime}(x_{1},x_{2})}\right]\cdot\mathbbm{1}{\left\\{\eta_{\Delta}(x_{1},x_{2})<0\right\\}}$
(50e) $\displaystyle=\begin{aligned}
&\Pr_{\bm{\epsilon}}\left[-\eta_{\Delta}(x_{1},x_{2})\cdot\epsilon^{\mathrm{AVE}}<\frac{\lvert
f(x_{1};\Delta)-f(x_{2};\Delta)\rvert}{K^{\frac{1}{\alpha}-1}\cdot\gamma^{\prime\prime}(x_{1},x_{2})}\right]\\\
&\cdot\mathbbm{1}{\left\\{\eta_{\Delta}(x_{1},x_{2})>0\right\\}}\\\
&+\Pr_{\bm{\epsilon}}\left[-\eta_{\Delta}(x_{1},x_{2})\cdot\epsilon^{\mathrm{AVE}}<\frac{\lvert
f(x_{1};\Delta)-f(x_{2};\Delta)\rvert}{K^{\frac{1}{\alpha}-1}\cdot\gamma^{\prime\prime}(x_{1},x_{2})}\right]\\\
&\cdot\mathbbm{1}{\left\\{\eta_{\Delta}(x_{1},x_{2})<0\right\\}}\end{aligned}$
(50f) $\displaystyle=\begin{aligned}
&\Pr_{\bm{\epsilon}}\left[-\eta_{\Delta}(x_{1},x_{2})\cdot\epsilon^{\mathrm{AVE}}<\frac{\lvert
f(x_{1};\Delta)-f(x_{2};\Delta)\rvert}{K^{\frac{1}{\alpha}-1}\cdot\gamma^{\prime\prime}(x_{1},x_{2})}\right]\\\
&\cdot\mathbbm{1}{\left\\{\eta_{\Delta}(x_{1},x_{2})\neq
0\right\\}}.\end{aligned}$ (50g)

This completes the proof. ∎

###  A-C Proof of Corollary 4

###### Proof.

Let $\beta_{m}=0,\forall m$, $\beta^{\prime\prime}(x_{1},x_{2})=0$, and
$\delta^{\prime\prime}(x_{1},x_{2})=0$ for any $x_{1},x_{2}$. We have

$-\eta_{\Delta}(x_{1},x_{2})\cdot\epsilon^{\mathrm{AVE}}\sim{\bf
S}(\alpha,0,1,0).$ (51)

This completes the proof. ∎

###  A-D Proof of Proposition 5

###### Proof.

Let $X_{1}=f(x_{1};{\bm{\epsilon}}_{1})-m_{f}(x_{1})$ and
$X_{2}=f(x_{2};{\bm{\epsilon}}_{2})-m_{f}(x_{2})$. They are independent, and
their distributions are symmetric around the origin. Let
$Z=\epsilon^{f}(x_{1},x_{2})-(m_{f}(x_{1})-m_{f}(x_{2}))$. Hence, we have
$Z=X_{1}-X_{2}$. If $Z$ is symmetric around the origin,
$\epsilon^{f}(x_{1},x_{2})$ is symmetric around $m_{f}(x_{1})-m_{f}(x_{2})$,
and the center of the median of $\epsilon^{f}(x_{1},x_{2})$ is
$m_{f}(x_{1})-m_{f}(x_{2})$. Therefore, it is sufficient to demonstrate that
$Z$ is symmetric around the origin. On the other hand, it is well-known that a
sum of independent random variables that are both symmetric around the origin
is symmetric around the origin [19]. Therefore, the center of the median of
$\epsilon^{f}(x_{1},x_{2})$ is $m_{f}(x_{1})-m_{f}(x_{2})$. This completes the
proof. ∎

###  A-E Proof of Proposition 6

###### Proof.

Let $\xi_{x}=f(x;\bm{\epsilon}_{x})-m_{f}(x)$ and
$\xi_{y}=f(y;\bm{\epsilon}_{y})-m_{f}(y)$, where $\bm{\epsilon}_{x}$ and
$\bm{\epsilon}_{y}$ are i.i.d. Then, $\xi_{x}$ and $\xi_{y}$ are i.i.d. as
well. Because they are i.i.d., $\xi_{x}-\xi_{y}$ is symmetric around $0$.
Therefore, the center of the median of $\xi_{x}-\xi_{y}$ is $0$. Now, consider
the center of the median of $f(x;\bm{\epsilon}_{x})-f(y;\bm{\epsilon}_{y})$,
namely $m_{f}(x,y)$. Because
$f(x;\bm{\epsilon}_{x})-f(y;\bm{\epsilon}_{y})=m_{f}(x)-m_{f}(y)+(\xi_{x}-\xi_{y})$,
$m_{f}(x,y)$ is the sum of $m_{f}(x)-m_{f}(y)$ and the center of the median of
$\xi_{x}-\xi_{y}$, the latter of which is $0$. Hence, we have
$m_{f}(x,y)=m_{f}(x)-m_{f}(y)$. ∎

###  A-F Proof of Proposition 7

###### Proof.

Let $X=f(x;\bm{\epsilon})-m_{f}(x)$. Then, $X$ is a random variable and the
center of its median is $0$. Because $\xi$ has a positive probability density
over $\mathbb{R}$, it has a unique median $m_{\xi}$. Let $Y=\xi-m_{\xi}$.
Obviously, its unique median is $0$. Because $X$ and $Y$ are independent and
$Y$ is absolutely continuous, $Z=X+Y$ is absolutely continuous. Let $F_{x}$ be
the cumulative distribution function of $X$ and $p_{y}$ be the probability
density function of $Y$. Let the probability density function of $Z$ be
denoted by $p_{z}$. Let $(L,U]$ be an interval with $F_{x}(U)-F_{x}(L)>0$.
Then, for any $z\in\mathbb{R}$, we have $\min_{t\in[z-U,z-L]}p_{y}(t)>0$ as
$p_{y}$ is continuous and strictly positive, and

$\displaystyle p_{z}(\xi)$
$\displaystyle=\int_{x\in\mathbb{R}}p_{y}(z-x)\mathrm{d}F_{x}(x)$ (52a)
$\displaystyle\geqslant\int_{x\in(L,U]}p_{y}(z-x)\mathrm{d}F_{x}(x)$ (52b)
$\displaystyle\geqslant\min_{t\in[z-U,z-L]}p_{y}(t)(F_{x}(U)-F_{x}(L))$ (52c)
$\displaystyle>0\enspace.$ (52d)

Therefore, $p_{z}(z)>0$ for all $z\in\mathbb{R}$. Because
$f(x;\epsilon)=Z+m_{f}(x)+m_{\xi}$, its probability density is
$p(t)=p_{z}(t-m_{f}(x)+m_{\xi})$. Then, it is easy to see that $p(t)>0$ for
all $t\in\mathbb{R}$. Because $p$ is the derivative of the cumulative
distribution function of $f(x;\epsilon)$, the cumulative distribution function
is differentiable. Therefore, 4 is satisfied. ∎

###  A-G Proof of Proposition 8

###### Proof.

We first show that (26) under 3. Let
$Z=\epsilon^{f}(x_{1},x_{2})-m_{f}(x_{1},x_{2})$. We have

$\displaystyle\mathbb{E}\left[\mathrm{sign}(\epsilon^{f}(x_{1},x_{2}))\right]$
$\displaystyle=\Pr\left[\epsilon^{f}(x_{1},x_{2})>0\right]-\Pr\left[\epsilon^{f}(x_{1},x_{2})<0\right]$
(53a)
$\displaystyle=\Pr\left[Z>-m_{f}(x_{1},x_{2})\right]-\Pr\left[Z<-m_{f}(x_{1},x_{2})\right].$
(53b)

Clearly, the first and second terms are no smaller and greater than
$\frac{1}{2}$, respectively, if $m_{f}(x_{1},x_{2})>0$. Under 3, we have
$m_{f}(x_{1},x_{2})>0\Leftrightarrow m_{f}(x_{1})>m_{f}(x_{2})$. Therefore,
$\mathbb{E}\left[\mathrm{sign}(\epsilon^{f}(x_{1},x_{2}))\right]\geqslant 0$
if $m_{f}(x_{1})>m_{f}(x_{2})$. Similarly, if $m_{f}(x_{1})<m_{f}(x_{2})$,
$\mathbb{E}\left[\mathrm{sign}(\epsilon^{f}(x_{1},x_{2}))\right]\leqslant 0$.
Thus, it is straightforward to obtain (26).

Next, we prove (27). Let $X_{1}=f(x_{1};{\bm{\epsilon}}_{1})-m_{f}(x_{1})$ and
$X_{2}=f(x_{2};{\bm{\epsilon}}_{2})-m_{f}(x_{2})$. Then, we have
$Z=X_{1}-X_{2}$. Let $\phi_{1}$, $\phi_{2}$, and $\phi_{z}$ be the
distribution functions of $X_{1}$, $X_{2}$, and $Z$, respectively. Suppose
that $\phi_{z}(z)>\phi_{z}(0)=\frac{1}{2}$ and
$\phi_{z}(-z)<\phi_{z}(0)=\frac{1}{2}$ for any $z>0$. Then, if
$m_{f}(x_{1},x_{2})>0$, we have

$\Pr\left[Z>-m_{f}(x_{1},x_{2})\right]=1-\phi_{z}(-m_{f}(x_{1},x_{2}))>\frac{1}{2}$
(54)

and

$-\Pr\left[Z<-m_{f}(x_{1},x_{2})\right]\leqslant\phi_{z}(-m_{f}(x_{1},x_{2}))\leqslant\frac{1}{2}.$
(55)

Thus, $\mathbb{E}\left[\mathrm{sign}(\epsilon^{f}(x_{1},x_{2}))\right]>0$.
Analogously, if $m_{f}(x_{1},x_{2})<0$,
$\mathbb{E}\left[\mathrm{sign}(\epsilon^{f}(x_{1},x_{2}))\right]<0$. Hence, it
suffices to show that $\phi_{z}(z)>\phi_{z}(0)=\frac{1}{2}$ and
$\phi_{z}(-z)<\phi_{z}(0)=\frac{1}{2}$ for any $z>0$ under 3 and 4. We already
know that $Z$ is symmetric around the origin and that the median is located at
the origin. Therefore, we have $\phi_{z}(z)\geqslant\frac{1}{2}$ for any
$z\geqslant 0$, $\phi(z)\leqslant\frac{1}{2}$ for $z<0$, and
$\phi_{z}(z)=1-\lim_{\rho\downarrow 0}\phi_{z}(-z-\rho)$. In other words, if
we show $\phi(z)<\frac{1}{2}$ for any $z<0$, then
$\phi(z)=1-\lim_{\rho\downarrow 0}\phi_{z}(-z-\rho)\geqslant
1-\phi(-z)>\frac{1}{2}$ for any $z>0$.

We prove that $\phi(z)<\frac{1}{2}$ for any $z<0$ holds. Under 4, there exists
$\delta>0$ such that $\phi_{1}$ and $\phi_{2}$ are continuous and strictly
increasing in interval $[-\delta,\delta]$. Let $-\delta\leqslant z<-\rho<0$.
We have

$\displaystyle\phi_{z}(z)$
$\displaystyle=\int^{\infty}_{-\infty}\phi_{1}(z-t)\mathrm{d}\phi_{2}(t)$
(56a) $\displaystyle=\begin{aligned}
&\int^{\infty}_{-\infty}\phi_{1}(-\rho-t)\mathrm{d}\phi_{2}(t)\\\
&+\int^{\infty}_{-\infty}\left(\phi_{1}(z-t)-\phi_{1}(-\rho-t)\right)\mathrm{d}\phi_{2}(t)\end{aligned}$
(56b)
$\displaystyle=\phi_{z}(-\rho)+\int^{\infty}_{-\infty}\left(\phi_{1}(z-t)-\phi_{1}(-\rho-t)\right)\mathrm{d}\phi_{2}(t)$
(56c)
$\displaystyle\leqslant\phi_{z}(-\rho)+\int_{z}^{-\rho}\left(\phi_{1}(z-t)-\phi_{1}(-\rho-t)\right)\mathrm{d}\phi_{2}(t)$
(56d)
$\displaystyle\leqslant\phi_{z}(-\rho)+\int_{z}^{-\rho}\left(\max_{z\leqslant
t\leqslant-\rho}\phi_{1}(z-t)-\phi_{1}(-\rho-t)\right)\mathrm{d}\phi_{2}(t)$
(56e) $\displaystyle=\begin{aligned} &\phi_{z}(-\rho)\\\
&+\left(\max_{z\leqslant
t\leqslant-\rho}\phi_{1}(z-t)-\phi_{1}(-\rho-t)\right)\left(\phi_{2}(-\rho)-\phi_{2}(z)\right).\end{aligned}$
(56f)

Here, we have the Lebesgue–Stieltjes integral. For the first inequality, we
consider that $\phi_{1}(z-t)-\phi_{1}(-\rho-t)\leqslant 0$ for all
$t\in\mathbb{R}$ because $\phi_{1}$ increases at least weakly. For the second
inequality, we consider that $\phi_{1}(z-t)-\phi_{1}(-\rho-t)<0$ for all
$t\in[z,-\rho]$ because $-\delta<z+\rho\leqslant
z-t<-\rho-t\leqslant-\rho-z<\delta$ and $\phi_{1}$ is strictly increasing in
$[-\delta,\delta]$. We also consider that $\phi_{1}(z-t)-\phi_{1}(-\rho-t)$ is
continuous on $t\ in[z,-\rho]$ and its maximum exists in $[z,-\rho]$.

The desired inequality is derived by evaluating each component on the
rightmost side of the inequality. For the first term,
$\phi_{z}(-\rho)\leqslant\frac{1}{2}$. For the component of $\phi_{2}$, we
have $\phi_{2}(-\rho)-\phi_{2}(z)>0$ because $\phi_{2}$ is strictly increasing
in $[-\delta,\delta]$. The strict negativity of component $\phi_{1}$ is
derived as follows. For $z<t<-\rho$, we have
$\phi_{1}(z-t)<\phi_{1}(0)=\frac{1}{2}$ and
$\phi_{1}(-\rho-t)>\phi_{1}(0)=\frac{1}{2}$. For $t=z$,
$\phi_{1}(z-t)<\phi_{1}(0)=\frac{1}{2}$ and
$\phi_{1}(-\rho-t)=\phi_{1}(-\rho-z)>\phi_{1}(0)=\frac{1}{2}$. For $t=-\rho$,
$\phi_{1}(z-t)=\phi_{1}(z+\rho)<\frac{1}{2}$ and
$\phi_{1}(-\rho-t)=\phi_{1}(0)=\frac{1}{2}$. Thus,
$\phi_{1}(z-t)-\phi_{1}(-\rho-t)<0$ for all $t\in[z,-\rho]$. Altogether, we
have $\phi_{z}(z)<\frac{1}{2}$ for $z\in[-\delta,0)$, which immediately leads
to $\phi_{z}(z)<\frac{1}{2}$ for all $z<0$ from the monotonicity of the
cumulative distribution function. ∎

###  A-H Proof of Theorem 9

###### Proof.

Let $(x_{1},x_{2})$ be a pair of points in $\mathbb{R}^{D}$ that satisfy
$\eta_{\mathrm{med}}(x_{1},x_{2})>0$. In addition, let
$\epsilon^{f}_{k}(x_{1},x_{2})$, where $k=1,\ldots,K$, be i.i.d. copies of
$\epsilon^{f}(x_{1},x_{2})$. We have

$\displaystyle\Pr_{\bm{\epsilon}}\left[\widehat{\eta}_{\bm{\epsilon}}^{\mathrm{SA}}(x_{1},x_{2})=\eta_{\mathrm{med}}(x_{1},x_{2})\right]$
$\displaystyle=\Pr_{\bm{\epsilon}}\left[\left(\frac{1}{K}\sum_{k}^{K}\mathrm{sign}(f(x_{1};{\bm{\epsilon}}_{1,k})-f(x_{2};{\bm{\epsilon}}_{2,k}))\right)\eta_{\mathrm{med}}(x_{1},x_{2})>0\right]$
(57a)
$\displaystyle=\Pr_{\bm{\epsilon}}\left[\frac{1}{K}\sum_{k}^{K}\mathrm{sign}(\epsilon^{f}_{k}(x_{1},x_{2}))>0\right]$
(57b)
$\displaystyle=1-\Pr_{\bm{\epsilon}}\left[\frac{1}{K}\sum_{k}^{K}\mathrm{sign}(\epsilon^{f}_{k}(x_{1},x_{2}))\leqslant
0\right].$ (57c)

Next, we evaluate the upper bound of
$\Pr_{\bm{\epsilon}}\left[\frac{1}{K}\sum_{k}^{K}\mathrm{sign}(\epsilon^{f}_{k}(x_{1},x_{2}))\leqslant
0\right]$. As
$\left|\mathrm{sign}\left(\epsilon^{f}(x_{1},x_{2})\right)\right|\leqslant 1$
surely, both
$\mathbb{E}\left[\mathrm{sign}\left(\epsilon^{f}(x_{1},x_{2})\right)\right]$
and
$\operatorname{Var}\left[\mathrm{sign}\left(\epsilon^{f}(x_{1},x_{2})\right)\right]$
are finite. Under 3 and condition $\eta_{\mathrm{med}}(x_{1},x_{2})>0$ and
considering Proposition 8, we have
$\mathbb{E}[\mathrm{sign}\left(\epsilon^{f}(x_{1},x_{2})\right)]\geqslant 0$.
Hence,

$\displaystyle\Pr_{\bm{\epsilon}}\left[\frac{1}{K}\sum_{k}^{K}\mathrm{sign}\left(\epsilon^{f}_{k}(x_{1},x_{2})\right)\leqslant
0\right]$
$\displaystyle=\Pr_{\bm{\epsilon}}\Bigg{[}K\cdot\mathbb{E}\left[\mathrm{sign}\left(\epsilon^{f}(x_{1},x_{2})\right)\right]-\sum_{k}^{K}\mathrm{sign}\left(\epsilon^{f}_{k}(x_{1},x_{2})\right)$
(58a) $\displaystyle\enspace\enspace\geqslant
K\cdot\mathbb{E}[\mathrm{sign}\left(\epsilon^{f}(x_{1},x_{2})\right)]\Bigg{]}.$
(58b)

Considering Hoeffding’s inequality, the righthand side of the above inequality
is bounded from above by

$\exp\left(-\frac{K\cdot\left(\mathbb{E}\left[\mathrm{sign}\left(\epsilon^{f}(x_{1},x_{2})\right)\right]\right)^{2}}{2}\right).$
(59)

Therefore, for $\eta(x_{1},x_{2})>0$, we have

$\Pr_{\bm{\epsilon}}\left[\widehat{\eta}_{\bm{\epsilon}}^{\mathrm{SA}}(x_{1},x_{2})=\eta_{\mathrm{med}}(x_{1},x_{2})\right]\\\
\geqslant
1-\exp\left(-\frac{K\cdot\left(\mathbb{E}\left[\mathrm{sign}\left(\epsilon^{f}(x_{1},x_{2})\right)\right]\right)^{2}}{2}\right).$
(60)

Similarly, assuming that $\eta(x_{1},x_{2})<0$, we obtain the same inequality
to that obtained for $\eta(x_{1},x_{2})>0$. This completes the proof. ∎

## References

  * [1] Improving evolution strategies through active covariance matrix adaptation.  In 2006 IEEE international conference on evolutionary computation, pages 2814–2821. IEEE, 2006.
  * [2] Diagonal acceleration for covariance matrix adaptation evolution strategies.  Evolutionary computation, 28(3):405–435, 2020.
  * [3] Alan Agresti.  Analysis of ordinal categorical data, volume 656.  John Wiley & Sons, 2010.
  * [4] Ali Ahrari, Saber Elsayed, Ruhul Sarker, Daryl Essam, and Carlos A Coello Coello.  Revisiting implicit and explicit averaging for noisy optimization.  IEEE Transactions on Evolutionary Computation, 2022.
  * [5] Youhei Akimoto, Anne Auger, and Nikolaus Hansen.  An ode method to prove the geometric convergence of adaptive stochastic algorithms.  Stochastic Processes and their Applications, 145:269–307, 2022\.
  * [6] Dirk V Arnold and H-G Beyer.  A general noise model and its effects on evolution strategy performance.  IEEE Transactions on Evolutionary Computation, 10(4):380–391, 2006\.
  * [7] Dirk V Arnold and Hans-Georg Beyer.  Local performance of the ($\mu$/$\mu$i, $\lambda$)-es in a noisy environment.  In Foundations of Genetic Algorithms 6, pages 127–141. Elsevier, 2001.
  * [8] Anne Auger and Nikolaus Hansen.  Linear convergence of comparison-based step-size adaptive randomized search via stability of markov chains.  SIAM Journal on Optimization, 26(3):1589–1624, 2016.
  * [9] H-G Beyer and Bernhard Sendhoff.  Functions with noise-induced multimodality: a test for evolutionary robust optimization-properties and performance analysis.  IEEE Transactions on Evolutionary Computation, 10(5):507–526, 2006\.
  * [10] Hans-Georg Beyer.  Toward a theory of evolution strategies: Some asymptotical results from the (1,+ $\lambda$)-theory.  Evolutionary computation, 1(2):165–188, 1993.
  * [11] Hans-Georg Beyer.  Evolutionary algorithms in noisy environments: Theoretical issues and guidelines for practice.  Computer methods in applied mechanics and engineering, 186(2-4):239–267, 2000.
  * [12] Hans-Georg Beyer and Bernhard Sendhoff.  Robust optimization–a comprehensive survey.  Computer methods in applied mechanics and engineering, 196(33-34):3190–3218, 2007.
  * [13] George Casella and Roger L Berger.  Statistical inference.  Cengage Learning, 2021.
  * [14] Benjamin Doerr and Andrew M Sutton.  When resampling to cope with noise, use median, not mean.  In Proceedings of the Genetic and Evolutionary Computation Conference, pages 242–248, 2019.
  * [15] Yinpeng Dong, Hang Su, Baoyuan Wu, Zhifeng Li, Wei Liu, Tong Zhang, and Jun Zhu.  Efficient decision-based black-box adversarial attacks on face recognition.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 7714–7722, 2019.
  * [16] Garuda Fujii, Masayuki Takahashi, and Youhei Akimoto.  Cma-es-based structural topology optimization using a level set boundary expression―application to optical and carpet cloaks.  Computer Methods in Applied Mechanics and Engineering, 332:624–643, 2018.
  * [17] Virginie Gabrel, Cécile Murat, and Aurélie Thiele.  Recent advances in robust optimization: An overview.  European journal of operational research, 235(3):471–483, 2014\.
  * [18] Saeed Ghadimi and Guanghui Lan.  Stochastic first-and zeroth-order methods for nonconvex stochastic programming.  SIAM Journal on Optimization, 23(4):2341–2368, 2013.
  * [19] Geoffrey Grimmett and David Stirzaker.  Probability and random processes.  Oxford university press, 2020.
  * [20] Nikolaus Hansen.  The cma evolution strategy: a comparing review.  Towards a new evolutionary computation: Advances in the estimation of distribution algorithms, pages 75–102, 2006.
  * [21] Nikolaus Hansen, Sibylle D Müller, and Petros Koumoutsakos.  Reducing the time complexity of the derandomized evolution strategy with covariance matrix adaptation (cma-es).  Evolutionary computation, 11(1):1–18, 2003.
  * [22] Nikolaus Hansen, André SP Niederberger, Lino Guzzella, and Petros Koumoutsakos.  A method for handling uncertainty in evolutionary optimization with an application to feedback control of combustion.  IEEE Transactions on Evolutionary Computation, 13(1):180–197, 2008\.
  * [23] Nikolaus Hansen and Andreas Ostermeier.  Completely derandomized self-adaptation in evolution strategies.  Evolutionary computation, 9(2):159–195, 2001.
  * [24] Yaochu Jin and Jürgen Branke.  Evolutionary optimization in uncertain environments-a survey.  IEEE Transactions on evolutionary computation, 9(3):303–317, 2005\.
  * [25] Iris Kriest, Volkmar Sauerland, Samar Khatiwala, Anand Srivastav, and Andreas Oschlies.  Calibrating a global three-dimensional biogeochemical ocean model (mops-1.0).  Geoscientific Model Development, 10(1):127–154, 2017.
  * [26] John A Nelder and Roger Mead.  A simplex method for function minimization.  The computer journal, 7(4):308–313, 1965.
  * [27] John P Nolan.  Univariate stable distributions.  Springer, 2020.
  * [28] Chao Qian, Yang Yu, Ke Tang, Yaochu Jin, Xin Yao, and Zhi-Hua Zhou.  On the effectiveness of sampling for evolutionary optimization in noisy environments.  Evolutionary computation, 26(2):237–267, 2018.
  * [29] Pratyusha Rakshit, Amit Konar, and Swagatam Das.  Noisy evolutionary optimization algorithms–a comprehensive survey.  Swarm and Evolutionary Computation, 33:18–45, 2017.

##  Appendix B Biography Section