# ProGO: Probabilistic Global Optimizer

Xinyu Zhang & Sujit Ghosh
Department of Statistics
NC State University
Raleigh, NC 27695, USA
<EMAIL_ADDRESS>

###### Abstract

In the field of global optimization, many existing algorithms face challenges
posed by non-convex target functions and high computational complexity or
unavailability of gradient information. These limitations, exacerbated by
sensitivity to initial conditions, often lead to suboptimal solutions or
failed convergence. This is true even for Metaheuristic algorithms designed to
amalgamate different optimization techniques to improve their efficiency and
robustness. To address these challenges, we develop a sequence of
multidimensional integration-based methods that we show to converge to the
global optima under some mild regularity conditions. Our probabilistic
approach does not require the use of gradients and is underpinned by a
mathematically rigorous convergence framework anchored in the nuanced
properties of nascent optima distribution. In order to alleviate the problem
of multidimensional integration, we develop a latent slice sampler that enjoys
a geometric rate of convergence in generating samples from the nascent optima
distribution, which is used to approximate the global optima. The proposed
Probabilistic Global Optimizer (ProGO) provides a scalable unified framework
to approximate the global optima of any continuous function defined on a
domain of arbitrary dimension. Empirical illustrations of ProGO across a
variety of popular non-convex test functions (having finite global optima)
reveal that the proposed algorithm outperforms, by order of magnitude, many
existing state-of-the-art methods, including gradient-based, zeroth-order
gradient-free, and some Bayesian Optimization methods, in term regret value
and speed of convergence. It is, however, to be noted that our approach may
not be suitable for functions that are expensive to compute.

##  1 Introcution

Global optimization constitutes a critical research area within applied
mathematics and numerical analysis, aiming to locate the global optima of
target functions over a specified domain. This field has substantial
applications across various sectors in machine learning, such as
hyperparameter tuning (Snoek et al., 2012), signal processing (Liu et al.,
2020), and black-box adversarial attacks (Ru et al., 2019). A global
optimization problem (minimization) with unique minima can be formulated as

${\bm{x}}^{*}=\operatorname*{arg\,min}_{{\bm{x}}\in\Omega}f({\bm{x}}),$ (1)

where $f({\bm{x}})$ is generally a continuous function defined over a domain
$\Omega$ $(\subseteq\mathbb{R}^{d})$, that is a subset of the $d$-dimensional
Euclidean space $\mathbb{R}^{d}$. Extensive progress has been made in
optimizing globally convex functions over compact domains, where the global
minima ${\bm{x}}^{*}$ is guaranteed to be identified. However, less
generalizable solutions exist for non-convex functions or on non-compact sets,
even when some target function possesses smoothness or differentiability.

For semantic precision, we differentiate between “optimum” / “minimum”, the
optimal / lowest function value $f^{*}$, and “optima” / “minima”, which
corresponds to the ${\bm{x}}^{*}$ at which $f^{*}$ is attained. In this paper,
we assume the existence of a finite
$f^{*}=\min_{{\bm{x}}\in\Omega}f({\bm{x}})$ and a non-empty set of minima
$\Omega^{*}=\\{{\bm{x}}\in\Omega:f({\bm{x}})=f^{*}\\}$. When the minima
${\bm{x}}^{*}$ in eq. (1) is unique, $\Omega^{*}$ will be a singleton set.
Importantly, $\Omega^{*}$ can comprise either a finite or infinite number of
elements. One of our main contributions lies in the identification of the
limitations of gradient-based techniques, particularly for non-convex
functions, and the introduction of a reliable, integration-based alternative
that guarantees to locate the global optima without convex assumptions.

Additionally, the efficacy of many global optimization algorithms is sensitive
to initial points. Even metaheuristic algorithms, which amalgamate various
optimization techniques for robustness, can yield suboptimal outcomes with
poorly chosen initial conditions. Our method, under mild regularity
conditions, is robust to initial conditions and yields accurate estimates of
$x^{*}$ within a decent computational timeframe, provided that function
evaluations are not expensive.

#### Gradient-based algorithms.

Gradient-based methods like stochastic gradient descent (Robbins & Monro,
1951), Adam (Kingma & Ba, 2015), AdaGrad (Duchi et al., 2011), and RMSprop
(Tieleman & Hinton, 2012) have found broad applications across disciplines.
They have demonstrated their utility in various successful applications such
as generative adversarial networks (Seward et al., 2018) and reinforcement
learning (Mnih et al., 2016). While these methods offer practical effective
solutions, their theoretical convergence to global optima is often framed
within specific contexts, particularly when the target function $f({\bm{x}})$
is smooth and convex. Recent variants like AMSGrad (Reddi et al., 2019) and
parameter-selective approaches (Shi et al., 2020) explore parameter effects on
theoretical guarantee and practical efficacy.

#### Zeroth order (ZO) methods.

In cases where gradient information is unavailable, noisy, or computationally
expensive to evaluate, such as signal processing and machine learning (Liu et
al., 2020), ZO methods have emerged driven by the need to solve these
problems. These techniques, also known as “black-box” or “derivative-free”
optimization, bypass the need for gradients and focus solely on function
values at any given point (Larson et al., 2019; Rios & Sahinidis, 2013).
Several noteworthy methods have been proposed in this vein, including
Gradientless Descent (GLD) (Golovin et al., 2019) which is numerically stable
via a geometric approach, Random Gradient-Free (RGF) method via finite
difference along a random direction by Nesterov & Spokoiny (2017), and Prior-
Guided Random Gradient-Free (PRGF) by Cheng et al. (2021), typically operating
under a framework that assumes convexity in the target functions.
Additionally, Shu et al. (2022) introduced the Zeroth-Order Optimization with
the Trajectory-Informed Derivative Estimation (ZoRD) algorithm, further
enriching the query-efficient ZO optimization methods landscape.

#### Global optimization challenges.

While gradient-based methods and ZO methods offer advantages, each comes with
its set of limitations and is often contingent upon wise selections of initial
parameters and starting points. Besides, Monte Carlo-based methods provide
consistent global convergence (Harrison, 2010) but can be computationally
demanding in high-dimensional spaces. Bayesian optimization techniques such as
the Gaussian Process-Upper Confidence Bound (GP-UCB) algorithm (Srinivas et
al., 2009) and Trust Region Bayesian Optimization (TuRBO) (Eriksson et al.,
2019) assume $f({\bm{x}})$ follows the Gaussian Process, whose performances
hinge upon careful selection of acquisition and kernel functions. Analytical
methods do contribute to domain-specific solutions but can entail intricate
numerical challenges, affecting their widespread applicability (Corriou &
Corriou, 2021). Most literature in this field has been oriented towards
establishing first-order optimality conditions, often under function convexity
and differentiability assumptions. A notable work by Luo (2018) formalized a
rigorous mathematical relation between an arbitrary continuous function $f$
defined over a compact set $\Omega\subseteq\mathbb{R}^{d}$ and its
corresponding global minima $f^{*}$; however, this work only built a
theoretical framework. This underscores the pressing need for an efficient and
robust global optimization framework, especially in addressing non-convex and
high-dimensional challenges.

(a) Surface of Ackley ($d=2$)  (b) Results on Ackley ($d=1000$)

Figure 1: (a) Visualization of the Ackley function for $d=2$. (b) Evaluation
of ProGO and competing methods (elaborated upon in Section 4) applied to the
Ackley function in a high-dimensional setting with $d=1000$. The $x$-axis
represents the iteration count, while the $y$-axis denotes the average log-
scaled regret. Each curve shows the mean $\pm$ standard error across ten
independent runs. ProGO outperforms other methods by a notably faster rate of
convergence accompanied by smaller variability.

#### Main contribution.

This paper introduces the Probabilistic Global Optimizer (ProGO), a novel non-
gradient-based global optimization algorithm based on a sequence of sampling
from a suitable probability distribution. Our work significantly extends the
theoretical framework laid by Luo (2018), notably in three key dimensions:

  1. 1.

Generalization to Non-Compact Set: Luo (2018)’s work is based on the
assumption that $\Omega$ is a compact set. Such an assumption may limit the
scope of its applicability to a class of popular functions when
$\Omega=\mathbb{R}^{d}$. E.g., even when $d=1$, the elementary function
$f({\bm{x}})={\bm{x}}^{2}$ defined over $\mathbb{R}$ is unbounded. We
generalize the domain to an arbitrary subset of $\mathbb{R}^{d}$ by
incorporating a sequence of probability distribution on $\mathbb{R}^{d}$ that
assigns less and less mass outside $\Omega^{*}$ and eventually converging with
support $\Omega^{*}$.

  2. 2.

A Practically Efficient Algorithm: While a theoretical framework is
established in Luo (2018) with several interesting results when the domain
$\Omega$ is compact and the function $f$ is assumed continuous or smooth (with
second order derivatives), to the best of our knowledge, no practical methods
are yet available to evaluate the (potentially high-dimensional) integration
required to estimate the minima. We fill this gap by developing a practically
efficient algorithm that we call ProGO, which uses a latent slice sampler
(explained later) to efficiently obtain samples from the probability
distribution of the minima.

  3. 3.

Extensive Experiments Validation: We carried out a comprehensive series of
experiments to evaluate ProGO’s performance in comparison to various types of
leading global optimization techniques, including Gradient Descent (GD),
Zeroth-Order Optimization (ZO), and Bayesian Optimization (BO). Our empirical
evidence demonstrates that ProGO consistently surpasses all the algorithms we
compared against across various metrics – most notably, geometric rate of
convergence to global minima and computational efficiency as indicated by
function evaluations or CPU time. Specifically, we illustrate the superior
numerical performance of our proposed ProGO for the popular Ackley function
(known to have several local optima with a global minimum at the origin) for
dimensions ranging from $d=20$ to $d=1000$ (refer to Fig. 1), as well as the
Levy function (refer to Fig. 5). As depicted in the figures, the logarithmic
regret exhibits a linear rate of convergence, which in the original scale
translates to geometric convergence, outpacing the majority of extant global
optimization algorithms.

The rest of this paper unfolds as follows: Section 2 lays the theoretical
groundwork on the probabilistic minima distribution for our approach; Section
3 details the latent slice sampler and our ProGO algorithm; Section 4 presents
the empirical validations, and Section 5 concludes the paper.

##  2 Preliminaries

###  2.1 Nascent Minima Probability Distribution

Our proposed algorithm, ProGO, is based on generalizing the sequence of
nascent minima (probability) distributions defined in Luo (2018) using an
arbitrary (prior) probability measure with full support on the Euclidean space
$\mathbb{R}^{d}$. This distribution possesses advantageous properties that
will be elaborated upon in subsequent discussions. In particular, we will show
how to efficiently generate samples from the sequence of such distributions
and subsequently use the empirical (posterior) mean and other summaries to
estimate the minimum value $f^{*}$ and the minima ${\bm{x}}^{*}$.

######  Assumption 1.

The following conditions are assumed throughout the paper:

  1. (i)

Assume that the function
$f:\Omega\subseteq\mathbb{R}^{d}\rightarrow\mathbb{R}$ is a continuous
function with a finite global minimum value $f^{*}$; i.e.,
$f({\bm{x}})\geqslant f^{*}$ for all ${\bm{x}}\in\Omega$.

  2. (ii)

The set of global minima
$\Omega^{*}=\\{{\bm{x}}\in\Omega:f({\bm{x}})=f^{*}\\}$ is non-empty.

  3. (iii)

There is a probability measure with density $\pi({\bm{x}})$ that has full
suppprt on $\mathbb{R}^{d}$. In other words, $\pi({\bm{x}})>0$ for any
${\bm{x}}\in\mathbb{R}^{d}$ and
$\int_{\mathbb{R}^{d}}\pi({\bm{x}})d{\bm{x}}=1$. Here, the integration is with
respect to the Lebesgue measure on $\mathbb{R}^{d}$.

In the above, the probability density $\pi({\bm{x}})$ can be chosen
arbitrarily, but in practice, we can use a uniform distribution when $\Omega$
is compact, or a very flat (nearly uniform) distribution when $\Omega$ is
unbounded. Regardless, we next define a nascent minima distribution when it
depends on the choice of the density $\pi(\cdot)$.

######  Definition 1.

Nascent Minima (probability) distribution:

For any $k\geqslant 0$, a nascent minima distribution density is defined as:

$\displaystyle
m_{k}({\bm{x}})=\frac{e^{-kf({\bm{x}})}\cdot\pi({\bm{x}})}{\int_{\Omega}e^{-kf({\bm{t}})}\cdot\pi({\bm{t}})\mathop{}\\!\mathrm{d}{\bm{t}}}.$
(2)

######  Remark 1.

Note that the denominator in eq. (2) is a finite positive quantity for any
arbitrary $k\geqslant 0$, because
$0<\int_{\Omega}e^{-kf({\bm{t}})}\cdot\pi(t)\mathop{}\\!\mathrm{d}t\leqslant
e^{-kf^{*}}$. The assumption that $\pi(\cdot)$ is a probability density can be
relaxed as long as
$\int_{\Omega}e^{-kf({\bm{t}})}\cdot\pi({\bm{t}})\mathop{}\\!\mathrm{d}{\bm{t}}<\infty$
for any $k>0$, even when $\int_{\Omega}\pi(t)\mathop{}\\!\mathrm{d}t=\infty$.

Notice that $m_{k}({\bm{x}})$ can also written as

$\displaystyle m_{k}({\bm{x}})$
$\displaystyle=\frac{e^{-k(f({\bm{x}})-f^{*})}\cdot\pi({\bm{x}})}{\int_{\Omega}e^{-k(f({\bm{t}})-f^{*})}\cdot\pi({\bm{t}})\mathop{}\\!\mathrm{d}{\bm{t}}}.$
(3)

By replacing the original $f({\bm{x}})$ by $f({\bm{x}})-f^{*}$, we may assume
without loss of generality that $f$ is a non-negative valued function and has
a global minimum $f^{*}=0$. Consequently, we focus on finding a solution in
the set of global minima
$\Omega^{*}=\\{{\bm{x}}\in\Omega:f({\bm{x}})=0\\}\neq\emptyset$ for the rest
of this paper for all subsequent theoretical analyses.

However, it should be noted that for practical applications, we need to work
with the original function $f$ as we will not know that global minimum $f^{*}$
and set of minima $\Omega^{*}$, and our goal would be to approximate $f^{*}$
and ${\bm{x}}^{*}\in\Omega^{*}$ by letting $k\rightarrow\infty$.

######  Remark 2.

If the original $f$ is a positive valued function with $f^{*}>0$, we can
replace it with $\log f$ when defining the nascent minima density in eq. (2).
Also, if a global maximum is desired, we can replace the original $f$ by $-f$
in defining the nascent minima density in eq. (2).

Next, we provide two results under very minimal conditions, which establish
the convergence of the (generalized) moments of the nascent minima
distribution to the minimum value $f^{*}$. With additional conditions, we also
establish the convergence of the minima.

###  2.2 Convergence Properties

In this subsection, we establish the theoretical underpinnings regarding the
convergence properties of the nascent minima distribution.

######  Theorem 1.

Consider a function $f$ and a probability density $\pi$ satisfying the
assumptions (i)-(iii) given in Assumption 1. Then, the nascent minima
distribution has the following properties:

$\lim_{k\rightarrow\infty}\int_{\Omega}f({\bm{x}})m_{k}({\bm{x}})\mathrm{d}{\bm{x}}=f^{*}$
(4)

The proof is detailed in Section A.1. The above result implies that the
expected value $\mathbb{E}_{m_{k}}[f(X)]$ converges to the minimum value
$f^{*}$ and hence, if we are able to generate samples from the nascent minima
distribution $m_{k}(\cdot)$ for any $k>0$, then we can approximate $f^{*}$
arbitrarily close by choosing a large $k>0$. Next, we show that the above
convergence is monotonic, which in turn implies that by increasing $k$
sequentially, we will get closer and closer to the minimum value $f^{*}$.

######  Theorem 2 (monotonicity).

Consider a non-constant function $f$ and a probability density $\pi$
satisfying the assumptions (i)-(iii) given in Assumption 1. For each $k>0$,
let $\mu_{k}=\mathbb{E}_{m_{k}}[f(X)]$ denotes the expectation of $f(X)$ when
$X\sim m_{k}(\cdot)$. Then the sequence $\\{\mu_{k}\\}$ is monotonically
decreasing and satisfies

$\displaystyle\frac{\mathrm{d}\mu_{k}}{\mathrm{d}k}=\frac{\mathrm{d}}{\mathrm{d}k}\int_{\Omega}f({\bm{x}})m_{k}({\bm{x}})\mathrm{d}{\bm{x}}=-\mathbb{V}{ar}_{k}(f)<0,$
(5)

where
$\mathbb{V}{ar}_{k}(f)=\int_{\Omega}\left(f({\bm{x}})-\mu_{k}\right)^{2}m_{k}({\bm{x}})\mathrm{d}{\bm{x}}$
denotes the variance of $f(X)$ when $X\sim m_{k}(\cdot)$.

The detailed proof is provided in Section A.2. Notice that the monotonic
convergence of $\mu_{k}$ as $k$ increases is established for any continuous
function $f$ and probability density $\pi$ satisfying (i)-(iii) of Assumption
1. This allows us for a very general use of ProGO with minimal assumptions for
any dimension $d\geqslant 1$. Next, we explore the convergence of the minimum
values with additional assumptions.

Here, we introduce the strong separability condition to define the scenario in
which this ProGo method is most suitable.

######  Assumption 2 (strong separability condition).

Consider the set of minima
$\Omega^{*}=\\{{\bm{x}}\in\Omega:f({\bm{x}})=f^{*}\\}$ which is assumed to be
non empty. Then $f$ is said to satisfy a strong separability if, for any given
$\delta>0$, we have
$\inf_{{\bm{x}}\notin\Omega^{*}:||{\bm{x}}-\tilde{{\bm{x}}}||>\delta}f({\bm{x}})>f^{*}$
for any $\tilde{{\bm{x}}}\in\Omega^{*}$.

The above condition implies that if ${\bm{x}}\notin\Omega^{*}$ and
$||{\bm{x}}-\tilde{{\bm{x}}}||>\delta$ for some $\delta>0$, then
$f({\bm{x}})>f^{*}$. In other words, the $f$ values for ${\bm{x}}$ not in
$\Omega^{*}$ are well separated from those that are in $\Omega^{*}$ and hence
$\epsilon_{0}=\inf\\{f({\bm{x}})-f^{*}:{\bm{x}}\notin\Omega^{*}\\}>0$.

######  Theorem 3.

Consider a bounded probability density $\pi$ and a target function satisfying
the assumptions (i)-(iii) in Assumption 1. For each $k>0$, let
$\Omega_{k}=\\{{\bm{x}}\in\Omega:m_{k}({\bm{x}})\geqslant
m_{k}(\tilde{{\bm{x}}}),\forall\tilde{{\bm{x}}}\in\Omega\\}$ denotes the set
of maximizes for $m_{k}({\bm{x}})$. Then for any sequence
$\\{{\bm{x}}_{k}^{*}\\}$, ${\bm{x}}_{k}^{*}\in\Omega_{k}$, it satisfy

$\lim_{k\rightarrow\infty}f({\bm{x}}_{k}^{*})=f^{*},$ (6)

In addition, if the target function satisfies the strong separability
condition in Assumption 2, we have

$\lim_{k\rightarrow\infty}\inf_{\tilde{{\bm{x}}}\in\Omega^{*}}||{\bm{x}}_{k}^{*}-\tilde{{\bm{x}}}||=0.$
(7)

The proof is available in Section A.3. This theorem establishes the
convergence of the sequence ${{\bm{x}}_{k}^{*}}$ toward the global minima set
$\Omega^{*}$ by quantifying the metric
$\min_{\tilde{{\bm{x}}}\in\Omega^{*}}||{\bm{x}}_{k}^{*}-\tilde{{\bm{x}}}||$.
This metric serves as a measure of divergence between the iteratively obtained
minima ${\bm{x}}_{k}^{*}$ for each $k$ and elements from the true global
minima set $\Omega^{*}$. The empirical evaluations in Section 4 corroborate
the algorithm’s efficacy in identifying discrepancies in both the optimal
value and corresponding minima.

##  3 ProGO: A New Probabilistic Global Optimization Method

This section describes Latent Slice Sampler (LSS) in Section 3.1, followed by
the ProGO algorithm in Section 3.2.

###  3.1 Latent Slice Sampler

Compared to traditional Markov Chain Monte Carlo techniques like the
Metropolis–Hastings algorithm, slice sampling (SS) provides benefits such as
reduced asymptotic variance and accelerated convergence (Mira & Tierney, 2002;
Neal, 2003; Roberts & Rosenthal, 1999). From several SS variants such as
Elliptical SS by Murray et al. (2010) and Polar SS by Rudolf & Schär (2023),
we adopt the LSS introduced by Li & Walker (2023) because it eliminates the
requirement for proposal distribution and improves efficiency in high-
dimensional sampling.

For a detailed understanding of LSS, consider the target distribution as a
minima distribution $m({\bm{x}})$ for a $d$-dimensional variable ${\bm{x}}$.
By incorporating slice variables $w$, ${\bm{s}}=(s_{1},\cdots,s_{d})^{\top}$,
and ${\bm{l}}=(l_{1},\cdots,l_{d})^{\top}$, the joint density can be
formulated as follow:

$p({\bm{x}},w,{\bm{s}},{\bm{l}})=\mathbb{I}\left(m({\bm{x}})>w\right)p({\bm{s}})\prod_{j=1}^{d}\frac{\mathbb{I}\left({\bm{x}}_{j}-{\bm{s}}_{j}/2<{\bm{l}}_{j}<{\bm{x}}_{j}+{\bm{s}}_{j}/2\right)}{{\bm{s}}_{j}}.$

Each ${\bm{s}}_{j}$ for $j=1,\cdots,d$ is assumed to follow an independent
gamma distribution with a shape parameter of 2 and scale parameter $\beta$ of
20, following Luo (2018). Let ${\bm{x}}^{(t)}$ represent the sample obtained
after $t^{th}$ iteration, the full LSS algorithm implemented via Gibbs
sampling is presented in Algorithm 1.

Input: Target probability distribution $\pi:\Omega\rightarrow[0,1]$, sample
size $N$, initialization ${\bm{x}}^{(0)}$, burn-in period $n_{b}$

1 Initialize ${\bm{x}}={\bm{x}}^{(0)}$, $w^{(0)}\sim
U(0,\pi({\bm{x}}^{(0)}))$, ${\bm{s}}^{(0)}_{j}\sim\mathrm{Gamma}(2,\beta)$,
$\mathrm{for}$ $j=1,\cdots,d,$ ${\bm{l}}^{(0)}\sim
U({\bm{x}}^{(0)}-{\bm{s}}^{(0)}/2,{\bm{x}}^{(0)}+{\bm{s}}^{(0)}/2)$

2 for  _$\mathrm{iteration}$ $t\in\\{1,2,\cdots,N+n_{b}\\}$_ do

3 ${\bm{a}}\leftarrow{\bm{l}}^{(t-1)}-{\bm{s}}^{(t-1)}/2$

4 ${\bm{b}}\leftarrow{\bm{l}}^{(t-1)}+{\bm{s}}^{(t-1)}/2$

5 while  _${\bm{x}}\notin\\{{\bm{x}}:\pi({\bm{x}}) >w^{(t-1)}\\}$_ do

6 for  _ $\mathrm{dimension}$ $j\in\\{1,\cdots,d\\}$_ do

7 if  _${\bm{x}} <{\bm{x}}_{j}$_ then

8 ${\bm{a}}_{j}\leftarrow\max({\bm{a}}_{j},{\bm{x}}_{j})$

9 else

10 ${\bm{b}}_{j}\leftarrow\min({\bm{b}}_{j},{\bm{x}}_{j})$

11 end if

12

13 end for

14 ${\bm{x}}\sim U({\bm{a}},{\bm{b}})$

15

16 end while

17 ${\bm{x}}^{(t)}\leftarrow{\bm{x}}$

18 $w^{(t)}\sim U(0,\pi({\bm{x}}^{(t)}))$

19 ${\bm{s}}^{(t)}\sim e^{\beta}+2|{\bm{l}}^{(t)}-{\bm{x}}^{(t)}|$

20 ${\bm{l}}^{(t)}\sim
U({\bm{x}}^{(t)}-{\bm{s}}^{(t)}/2,{\bm{x}}^{(t)}+{\bm{s}}^{(t)}/2)$

21

22 end for

23

Output: Samples $\\{{\bm{x}}^{(n_{b}+1)},\cdots,{\bm{x}}^{(n_{b}+N)}\\}$ from
target probability distribution $\pi$

Algorithm 1 Latent Slice Sampler for ProGO (LSS-ProGO)

#### One Dimensional Illustration.

To demonstrate the rationale behind using the nascent minima distribution for
global optimization and to evaluate the efficacy of LSS, we consider a one-
dimensional example adapted from Luo (2018), given by
$f(x)=\cos(x^{2})+x/5+1,x\in[0,5].$ This function possesses three local minima
and a singular global minima at $f(x)=0.353$ when $x=1.756$.

(a) $k=0$ (b) $k=1$

(c) $k=3$ (d) $k=9$

Figure 2: One dimensional illustration of nascent minima distribution
$m_{k}({\bm{x}})$ and its latent slice samples. In this illustration,
$\pi({\bm{x}})$ is assumed to satisfy Assumption 1 and uniformly distributed
within [0,5]. The target function $f(x)=\cos(x^{2})+x/5+1$ is depicted in the
blue curve. The true $m_{k}({\bm{x}})$ with $k=0,1,3,9$ are shown as a red
curve. The black histograms represent the latent slice samples of
$m_{k}({\bm{x}})$.

As depicted in Figure 2, the density distributions generated through LSS
closely align with the actual minima distribution across various $k$. Notably,
when $k=0$, $m_{0}({\bm{x}})=\pi({\bm{x}})$ corresponds to a uniform
distribution, as shown in Figure 2(a). A higher value of $k$ leads to a minima
distribution where the probability is increasingly focused on the global
optima.

###  3.2 ProGO

We outline the algorithm of ProGO in Algorithm 2. The step size is incremented
by $\triangle k$ for each iteration, chosen as $\triangle k=(e-1)k$. Such
choice is motivated by the inverse shrinkage rate of $m_{k}({\bm{x}})$ with
respect to $k$, as demonstrated in Luo (2018). The starting value of $k$
establishes the subsequent sequence $\\{k,ke,ke^{2},\cdots,ke^{T}\\}$ across
$T$ iterations, leading to a converging sequence of
$\\{f({\bm{x}}_{k}^{*}),f({\bm{x}}_{ke}^{*}),f({\bm{x}}_{ke^{2}}^{*}),\cdots,f({\bm{x}}_{ke^{T}}^{*})\\}$,
where $x_{k}^{*}\in\\{{\bm{x}}\in\Omega:m_{k}({\bm{x}})\geqslant
m_{k}(\tilde{{\bm{x}}}),\forall\tilde{{\bm{x}}}\in\Omega\\}$. Based on our
preliminary results, the initial value of $k$ is set to 5.

Input: Target function $f:\Omega\rightarrow\mathbb{R}$, probability
distribution function $\pi(\cdot)$, max iteration number $T=200$, sample size
$N$, starting point ${\bm{x}}^{(0)}$, burn-in period $n_{b}$

1 Initialize $k=5$

2 for  _$\mathrm{iteration}$ $t\in\\{1,2,\cdots,T\\}$_ do

3 ${\bm{x}}^{(1)},\cdots,{\bm{x}}^{(N)}\sim\text{LSS-
ProGO}(m_{k}(\cdot~{};\pi,f),N,{\bm{x}}^{(0)},n_{b})$

4
$\tilde{{\bm{x}}}^{(t)}\leftarrow\operatorname*{arg\,max}_{i\in\\{1,2,\cdots,N\\}}m_{k}({\bm{x}}^{(i)};\pi,f))$

5 $k\leftarrow k+(e-1)k$

6 end for

Output:
$\operatorname*{arg\,min}_{t\in\\{1,2,\cdots,T\\}}f(\tilde{{\bm{x}}}^{(t)})$

Algorithm 2 ProGO Algorithm

It is noteworthy that the output of ProGO delivers more than just an optimum
value; it also provides the sample sets from the minima distribution. This
provides valuable information about the distributional properties of local
minima, as exemplified in Figure 2.

##  4 Experiments

Test functions, evaluation metrics, and benchmarking methods play crucial
roles in validating the performance of an optimization algorithm. While a wide
array of test functions exists in the literature (Jamil & Yang, 2013), the
Ackley function—formulated initially by Ackley (2012)—and the Levy function
remain prevalent choices, as corroborated by recent studies like Shu et al.
(2022).

As for evaluation criteria, we use the following metrics of function log
regret and minima log regret to capture discrepancies in both $f({\bm{x}})$
and ${\bm{x}}$.

######  Definition 2.

Given $\tilde{{\bm{x}}}$ as an estimated optima in a $d$-dimension space, the
function log regret is defined as:

$r_{f}=\log\left(f(\tilde{{\bm{x}}})-f^{*}\right),$ (8)

quantifying the deviation between the estimated and true global optimum. The
minima log regret is formulated as:

$r_{m}=\log\frac{||\tilde{{\bm{x}}}-{\bm{x}}^{*}||}{\sqrt{d}},$ (9)

which quantifies the discrepancy between the estimated and true global minima.

The experimental design and the selection of competing algorithms are mostly
aligned with the framework presented in ZoRD (Shu et al., 2022) to ensure a
consistent and fair evaluation (details in Section B). A computational budget
capped at 200 iterations is allotted for each of the ten independent runs
conducted in every experiment setting. The evaluated methods include: 1) ZoRD:
Zeroth-order trajectory-informed derivative estimation (Shu et al., 2022). 2)
GD: Gradient-Descent, directly using first-order information. 3) GLD:
Gradientless Descent (Golovin et al., 2019). 4) PRGF: Prior-guided random
gradient-free algorithm (Cheng et al., 2021). 5) RGF: Random gradient-free
method via finite difference (Nesterov & Spokoiny, 2017). 6) GP-UCB: Gaussian
Process-Upper Confidence Bound (Srinivas et al., 2009). 7) TuRBO: Trust Region
Bayesian Optimization (Eriksson et al., 2019). 8) ProGO (our approach):
Probabilistic Global Optimization.

###  4.1 Ackley

The Ackley function serves as a prominent benchmark for evaluating
optimization algorithms. It is characterized as a continuous, differentiable,
multimodal, and non-convex function, thereby posing significant optimization
challenges. Mathematically, the $d$-dimensional Ackley function is given by:

$f({{\bm{x}}})=-a\exp\left(-b\sqrt{\frac{1}{d}\sum_{i=1}^{d}{\bm{x}}_{i}^{2}}\right)-\exp\left(\frac{1}{d}\sum_{i=1}^{d}\cos\left(c{\bm{x}}_{i}\right)\right)+a+\exp(1),$
(10)

where recommended parameter values are $a=20,b=0.2$ and $c=2\pi$ (Adorio &
Diliman, 2005). As shown in Figure 1 (a), the Ackley function features
numerous local minima and a central global minima at $f^{*}=0$ and
${\bm{x}}^{*}=(0,\cdots,0)^{\top}$, presenting multiple local minima traps for
optimization algorithms.

(a) Ackley (d=20) (b) Ackley (d=40)

Figure 3: Evaluation of ProGO and competing methods applied to the Ackley
function with $d=20,40$. The $x$-axis represents the iteration count, while
the $y$-axis denotes the average log-scaled regret. Each curve shows the mean
$\pm$ standard error across ten independent runs.  Table 1: Comparison of
performance on the Ackley function for dimensions $d=20,40,1000$: Each
optimization method has ten independent runs. Accuracy is quantified using the
average function log regret ($r_{f}$) and the average minima log regret
($r_{m}$) across these ten runs. Computational efficiency is represented by
the average runtime ($t$) across all ten runs, measured in seconds.

| $d=20$ |  | $d=40$ |  | $d=1000$
---|---|---|---|---|---
Method | ${r_{f}}$ | ${r_{m}}$ | ${t~{}(\mathrm{s})}$ |  | ${r_{f}}$ | ${r_{m}}$ | ${t~{}(\mathrm{s})}$ |  | ${r_{f}}$ | ${r_{m}}$ | ${t~{}(\mathrm{s})}$
ZoRD | 2.83 | 1.68 | 1265.43 |  | 2.85 | 1.72 | 1092.04 |  | 2.73 | 1.58 | 406.16
GD | 2.85 | 1.72 | 0.12 |  | 2.87 | 1.75 | 0.13 |  | 2.89 | 1.75 | 0.14
GLD | 2.73 | 1.65 | 0.38 |  | 2.89 | 1.70 | 0.39 |  | 2.97 | 1.75 | 0.40
PRGF | 2.85 | 1.70 | 0.06 |  | 2.89 | 1.74 | 0.07 |  | 2.97 | 1.75 | 0.08
RGF | 2.86 | 1.66 | 0.06 |  | 2.89 | 1.72 | 0.07 |  | 2.97 | 1.74 | 0.08
GP-UCB | 2.01 | 1.62 | 132.13 |  | 2.07 | 1.63 | 323.12 |  | 2.93 | 1.74 | 1132.61
TuRBO | 1.73 | 1.60 | 30.83 |  | 2.54 | 1.62 | 83.84 |  | 2.95 | 1.73 | 272.56
ProGO | -35.35 | -35.86 | 6.77 |  | -31.56 | -9.50 | 16.25 |  | 1.99 | 0.52 | 280.05

Our empirical evaluations (see Table 1) span dimensions $d=20,40,$ and $1000$.
Across all dimensions, ProGO consistently outperforms other methods, achieving
significantly lower function log regret and lower minima log regret. Moreover,
Figure 1(b) corroborates ProGO’s geometric rate of convergence even on high
dimensions (see Figure 3 for results on $d=20$ and $d=40$).

(a) Levy (d=40) (b) Levy (d=100)

Figure 4: Evaluation of ProGO and competing methods applied to the Ackley
function with $d=40,100$. The $x$-axis represents the iteration count, while
the $y$-axis denotes the average log-scaled regret. Each curve shows the mean
$\pm$ standard error across ten independent runs.

###  4.2 Levy

(a) Surface of Levy ($d=2$) (b) Results on Levy ($d=1000$)

Figure 5: (a) Visualization of the Levy function for $d=2$. (b) Evaluation of
ProGO and competing methods applied to the Levy function in a high-dimensional
setting with $d=1000$. The $x$-axis represents the iteration count, while the
$y$-axis denotes the average log-scaled regret. Each curve shows the mean
$\pm$ standard error across ten independent runs.

The Levy function is another frequently used test function in optimization
research, as in the work of Shu et al. (2022). It is a continuous and non-
convex function defined as:

$f({\bm{x}})=\sin^{2}\left(\pi
w_{1}\right)+\sum_{i=1}^{d-1}\left(w_{i}-1\right)^{2}\left[1+10\sin^{2}\left(\pi
w_{i}+1\right)\right]+\left(w_{d}-1\right)^{2}\left[1+\sin^{2}\left(2\pi
w_{d}\right)\right],$ (11)

where $w_{i}=1+({\bm{x}}_{i}-1)/4,\forall i=1,\ldots,d$. The global minimum is
$f^{*}=0$ attained at ${{\bm{x}}}^{*}=(1,\ldots,1)^{\top}$. As shown in Figure
5, the Levy function presents a more complex optimization challenge than the
Ackley function due to the substantially flatter area that surrounds its
global optima.

Table 2: Comparison of performance on the Levy function for dimensions
$d=40,100,1000$: Each optimization method undergoes ten independent runs.
Accuracy is quantified using the average function log regret ($r_{f}$) and the
average minima log regret ($r_{m}$) across these ten runs. Computational
efficiency is represented by the average runtime ($t$) across all ten runs,
measured in seconds.

| $d=40$ |  | $d=100$ |  | $d=1000$
---|---|---|---|---|---
Method | ${r_{f}}$ | ${r_{m}}$ | ${t~{}(\mathrm{s})}$ |  | ${r_{f}}$ | ${r_{m}}$ | ${t~{}(\mathrm{s})}$ |  | ${r_{f}}$ | ${r_{m}}$ | ${t~{}(\mathrm{s})}$
ZoRD | 4.08 | 1.56 | 585.69 |  | 5.26 | 1.56 | 1514.97 |  | 8.05 | 1.58 | 272.29
GD | 4.10 | 1.58 | 0.16 |  | 5.15 | 1.61 | 0.15 |  | 7.59 | 1.59 | 0.14
GLD | 4.89 | 1.59 | 0.40 |  | 6.15 | 1.60 | 0.39 |  | 8.88 | 1.59 | 0.39
PRGF | 4.68 | 1.59 | 0.08 |  | 5.96 | 1.61 | 0.07 |  | 8.75 | 1.60 | 0.07
RGF | 4.59 | 1.58 | 0.08 |  | 5.84 | 1.62 | 0.08 |  | 8.74 | 1.59 | 0.09
GP-UCB | 3.80 | 1.58 | 216.83 |  | 5.28 | 1.58 | 651.38 |  | 8.49 | 1.58 | 1400.41
TuRBO | 3.58 | 1.56 | 77.66 |  | 5.18 | 1.55 | 197.78 |  | 8.70 | 1.59 | 265.71
ProGO | -0.05 | -0.96 | 11.85 |  | 1.55 | -0.50 | 27.33 |  | 4.62 | -0.01 | 311.99

Empirical evaluation of the Levy function for dimensions $d=40,100,1000$ is
presented in Table 2. Notably, ProGO demonstrates significantly lower log
regrets relative to other competing methods across all dimensions.
Furthermore, Figure 5 and Figure 4 show that ProGO exhibits a markedly faster
convergence rate compared to competing gradient-based, zeroth-order, and
Bayesian optimization methods. Notice that Gradient Descent, depicted in
orange, initially exhibits rapid convergence but is then trapped in local
optima.

##  5 Conclusion

In this paper, we introduce ProGO, a novel probabilistic global optimization
algorithm leveraging minima distribution theory and latent slice sampling
technique. Our methodology represents a significant departure from
conventional gradient-based methods, offering robust convergence guarantees
for global optima while preserving computational efficiency without using
gradient information.

Specifically, our contributions are threefold: We extend Luo (2018)’s
theoretical framework to non-compact sets and prove its global convergence.
Based on the generalized framework, we implement the ProGO algorithm,
integrating a latent slice sampler for enhanced computational efficiency,
especially for high dimensions. Finally, comprehensive experiments demonstrate
ProGO’s outstanding performance over state-of-the-art methods in terms of
accuracy and convergence speed on various functions and dimensions.

However, it is worth noting that ProGO may not be suitable for optimization
problems where function evaluation is computationally expensive. Future
investigations on enhancing the algorithm’s computational efficiency and
extending the applicability of ProGO to more diverse problem domains may
contribute to the growing field of global optimization.

## References

  * (1)
  * Ackley (2012) Ackley, D. (2012), A connectionist machine for genetic hillclimbing, Vol. 28, Springer science & business media.
  * Adorio & Diliman (2005) Adorio, E. P. & Diliman, U. (2005), ‘Mvf-multivariate test functions library in c for unconstrained global optimization’, Quezon City, Metro Manila, Philippines 44.
  * Cheng et al. (2021) Cheng, S., Wu, G. & Zhu, J. (2021), ‘On the convergence of prior-guided zeroth-order optimization algorithms’, Advances in Neural Information Processing Systems 34, 14620–14631.
  * Corriou & Corriou (2021) Corriou, J.-P. & Corriou, J.-P. (2021), ‘Analytical methods for optimization’, Numerical methods and optimization: Theory and practice for engineers pp. 455–503.
  * Duchi et al. (2011) Duchi, J., Hazan, E. & Singer, Y. (2011), ‘Adaptive subgradient methods for online learning and stochastic optimization.’, Journal of machine learning research 12(7).
  * Eriksson et al. (2019) Eriksson, D., Pearce, M., Gardner, J., Turner, R. D. & Poloczek, M. (2019), ‘Scalable global optimization via local bayesian optimization’, Advances in neural information processing systems 32.
  * Golovin et al. (2019) Golovin, D., Karro, J., Kochanski, G., Lee, C., Song, X. & Zhang, Q. (2019), ‘Gradientless descent: High-dimensional zeroth-order optimization’, arXiv preprint arXiv:1911.06317 .
  * Harrison (2010) Harrison, R. L. (2010), Introduction to monte carlo simulation, in ‘AIP conference proceedings’, Vol. 1204, American Institute of Physics, pp. 17–21.
  * Jamil & Yang (2013) Jamil, M. & Yang, X.-S. (2013), ‘A literature survey of benchmark functions for global optimisation problems’, International Journal of Mathematical Modelling and Numerical Optimisation 4(2), 150–194.
  * Kingma & Ba (2015) Kingma, D. & Ba, J. (2015), Adam: A method for stochastic optimization, in ‘International Conference on Learning Representations (ICLR)’, San Diega, CA, USA.
  * Larson et al. (2019) Larson, J., Menickelly, M. & Wild, S. M. (2019), ‘Derivative-free optimization methods’, Acta Numerica 28, 287–404.
  * Li & Walker (2023) Li, Y. & Walker, S. G. (2023), ‘A latent slice sampling algorithm’, Computational Statistics & Data Analysis 179, 107652.
  * Liu et al. (2020) Liu, S., Chen, P.-Y., Kailkhura, B., Zhang, G., Hero III, A. O. & Varshney, P. K. (2020), ‘A primer on zeroth-order optimization in signal processing and machine learning: Principals, recent advances, and applications’, IEEE Signal Processing Magazine 37(5), 43–54.
  * Luo (2018) Luo, X. (2018), ‘Minima distribution for global optimization’, arXiv preprint arXiv:1812.03457 .
  * Mira & Tierney (2002) Mira, A. & Tierney, L. (2002), ‘Efficiency and convergence properties of slice samplers’, Scandinavian Journal of Statistics 29(1), 1–12.
  * Mnih et al. (2016) Mnih, V., Badia, A. P., Mirza, M., Graves, A., Lillicrap, T., Harley, T., Silver, D. & Kavukcuoglu, K. (2016), Asynchronous methods for deep reinforcement learning, in ‘International conference on machine learning’, PMLR, pp. 1928–1937.
  * Murray et al. (2010) Murray, I., Adams, R. & MacKay, D. (2010), Elliptical slice sampling, in ‘Proceedings of the thirteenth international conference on artificial intelligence and statistics’, JMLR Workshop and Conference Proceedings, pp. 541–548.
  * Neal (2003) Neal, R. M. (2003), ‘Slice sampling’, The annals of statistics 31(3), 705–767.
  * Nesterov & Spokoiny (2017) Nesterov, Y. & Spokoiny, V. (2017), ‘Random gradient-free minimization of convex functions’, Foundations of Computational Mathematics 17, 527–566.
  * Reddi et al. (2019) Reddi, S. J., Kale, S. & Kumar, S. (2019), ‘On the convergence of adam and beyond’, arXiv preprint arXiv:1904.09237 .
  * Rios & Sahinidis (2013) Rios, L. M. & Sahinidis, N. V. (2013), ‘Derivative-free optimization: a review of algorithms and comparison of software implementations’, Journal of Global Optimization 56, 1247–1293.
  * Robbins & Monro (1951) Robbins, H. & Monro, S. (1951), ‘A stochastic approximation method’, The annals of mathematical statistics pp. 400–407.
  * Roberts & Rosenthal (1999) Roberts, G. O. & Rosenthal, J. S. (1999), ‘Convergence of slice sampler markov chains’, Journal of the Royal Statistical Society Series B: Statistical Methodology 61(3), 643–660.
  * Ru et al. (2019) Ru, B., Cobb, A., Blaas, A. & Gal, Y. (2019), Bayesopt adversarial attack, in ‘International Conference on Learning Representations’.
  * Rudolf & Schär (2023) Rudolf, D. & Schär, P. (2023), ‘Dimension-independent spectral gap of polar slice sampling’, arXiv preprint arXiv:2305.03685 .
  * Seward et al. (2018) Seward, C., Unterthiner, T., Bergmann, U., Jetchev, N. & Hochreiter, S. (2018), First order generative adversarial networks, in ‘International Conference on Machine Learning’, PMLR, pp. 4567–4576.
  * Shi et al. (2020) Shi, N., Li, D., Hong, M. & Sun, R. (2020), Rmsprop converges with proper hyper-parameter, in ‘International Conference on Learning Representations’.
  * Shu et al. (2022) Shu, Y., Dai, Z., Sng, W., Verma, A., Jaillet, P. & Low, B. K. H. (2022), Zeroth-order optimization with trajectory-informed derivative estimation, in ‘The Eleventh International Conference on Learning Representations’.
  * Snoek et al. (2012) Snoek, J., Larochelle, H. & Adams, R. P. (2012), ‘Practical bayesian optimization of machine learning algorithms’, Advances in neural information processing systems 25.
  * Srinivas et al. (2009) Srinivas, N., Krause, A., Kakade, S. M. & Seeger, M. (2009), ‘Gaussian process optimization in the bandit setting: No regret and experimental design’, arXiv preprint arXiv:0912.3995 .
  * Tieleman & Hinton (2012) Tieleman, T. & Hinton, G. (2012), ‘Rmsprop: Divide the gradient by a running average of its recent magnitude. coursera: Neural networks for machine learning’, COURSERA Neural Networks Mach. Learn 17.

##  Appendix A Proofs

###  A.1 Proof for Theorem 1

###### Proof for Theorem 1.

Recall that $\Omega^{*}=\\{{\bm{x}}\in\Omega:f({\bm{x}})=0\\}$ without loss of
generalizability. For any $\epsilon>0$, define
$\Omega_{\epsilon}=\\{{\bm{x}}\in\Omega:f({\bm{x}})<\epsilon\\}$ and
$\Omega_{\epsilon}^{c}=\Omega\setminus\Omega_{\epsilon}$, then

$\displaystyle\int_{\Omega}f({\bm{x}})m_{k}({\bm{x}})d{\bm{x}}$
$\displaystyle=\frac{\int_{\Omega_{\frac{\epsilon}{2}}}f({\bm{x}})\cdot
e^{-kf({\bm{x}})}\pi({\bm{x}})\mathrm{d}{\bm{x}}+\int_{\Omega_{\frac{\epsilon}{2}}^{c}}f({\bm{x}})\cdot
e^{-kf({\bm{x}})}\cdot\pi({\bm{x}})\mathrm{d}{\bm{x}}}{\int_{\Omega}e^{-kf(t)}\cdot\pi(t)\mathrm{d}t}$
$\displaystyle<\frac{\epsilon}{2}\cdot\frac{\int_{\Omega_{\frac{\epsilon}{2}}}e^{-kf({\bm{x}})}\cdot\pi({\bm{x}})\mathrm{d}{\bm{x}}}{\int_{\Omega}e^{-kf(t)}\cdot\pi(t)\mathrm{d}t}+\frac{\int_{\Omega_{\frac{\epsilon}{2}}^{c}}f({\bm{x}})\cdot
e^{-kf({\bm{x}})}\cdot\pi({\bm{x}})\mathrm{d}{\bm{x}}}{\int_{\Omega_{\frac{\epsilon}{4}}}e^{-kf(t)}\cdot\pi(t)\mathrm{d}t}$
$\displaystyle\leqslant\frac{\epsilon}{2}\cdot
1+\frac{\int_{\Omega_{\frac{\epsilon}{2}}^{c}}f({\bm{x}})e^{-k\left(f({\bm{x}})-\frac{\epsilon}{4}\right)}\cdot\pi({\bm{x}})\mathrm{d}{\bm{x}}}{\int_{\Omega_{\frac{\epsilon}{4}}}\pi(t)\mathrm{d}t}$
$\displaystyle=\frac{\epsilon}{2}+\int_{\Omega_{\frac{\epsilon}{2}}^{c}}\frac{f({\bm{x}})\pi({\bm{x}})e^{-k\left(f({\bm{x}})-\frac{\epsilon}{4}\right)}}{\int_{\Omega_{\frac{\epsilon}{4}}}\pi(t)\mathrm{d}t}\mathrm{d}{\bm{x}}\stackrel{{\scriptstyle\Delta}}{{=}}\frac{\epsilon}{2}+\int_{\Omega_{\frac{\epsilon}{2}}^{c}}g_{k}({\bm{x}})\mathrm{d}{\bm{x}}.$

The first inequality lies in $f({\bm{x}})<\frac{\epsilon}{2}$ for any
${\bm{x}}\in\Omega_{\frac{\epsilon}{2}}$ and
$1/{\int_{\Omega}e^{-kf(t)}\cdot\pi(t)\mathrm{d}t}<1/{\int_{\Omega_{\frac{\epsilon}{2}}}e^{-kf(t)}\cdot\pi(t)\mathrm{d}t}$.
The second inequality is given by
${\int_{\Omega_{\frac{\epsilon}{2}}}e^{-kf({\bm{x}})}\cdot\pi({\bm{x}})\mathrm{d}{\bm{x}}}<{\int_{\Omega}e^{-kf({\bm{x}})}\cdot\pi({\bm{x}})\mathrm{d}{\bm{x}}}$
and $f({\bm{x}})<\frac{\epsilon}{4}$ for any
${\bm{x}}\in\Omega_{\frac{\epsilon}{4}}$. Define
$g_{k}({\bm{x}})=\frac{f({\bm{x}})\pi({\bm{x}})e^{-k\left(f({\bm{x}})-\frac{\epsilon}{4}\right)}}{\int_{\Omega_{\frac{\epsilon}{4}}}\pi(t)\mathrm{d}t}$
for ${\bm{x}}\in{\Omega_{\frac{\epsilon}{2}}^{c}}$, then
$\\{g_{k}({\bm{x}})\\}_{k=1}^{\infty}$ is a sequence of nonnegative functions
that monotonously decreases to 0 when k goes to infinity and converges to
zero, i.e., $\lim_{k\rightarrow\infty}g_{k}({\bm{x}})=0$. Hence, by monotone
convergence theorem,
$\lim_{k\rightarrow\infty}{\int_{\Omega_{\frac{\epsilon}{2}}^{c}}g_{k}({\bm{x}})\mathrm{d}{\bm{x}}}=0.$
Consequently, for any $\epsilon>0$, there exists a large $K_{0}$, such that
$\int_{\Omega_{\frac{\epsilon}{2}}^{c}}g_{k}({\bm{x}})\mathrm{d}{\bm{x}}<\epsilon/2$
for any $k>K_{0}$, and thus

$\displaystyle 0\leqslant\int_{\Omega}f({\bm{x}})m_{k}({\bm{x}})d{\bm{x}}$
$\displaystyle<\frac{\epsilon}{2}+\frac{\epsilon}{2}=\epsilon,$

which proves
$\lim_{k\rightarrow\infty}\int_{\Omega}f({\bm{x}})m_{k}({\bm{x}})\mathrm{d}{\bm{x}}=f^{*}(=0).$

∎

###  A.2 Proof for Theorem 2

###### Proof for Theorem 2.

For every $k\in\mathbb{R}$,

$\displaystyle\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}k}\log
m_{k}({\bm{x}})$
$\displaystyle=\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}k}\left(-kf({\bm{x}})+\log\pi({\bm{x}})-\log\left\\{\int_{\Omega}e^{-kf({\bm{t}})}\cdot\pi({\bm{t}})\mathop{}\\!\mathrm{d}{\bm{t}}\right\\}\right)$
$\displaystyle\stackrel{{\scriptstyle}}{{=}}-f({\bm{x}})-\frac{\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}k}\left\\{\int_{\Omega}e^{-kf({\bm{t}})}\cdot\pi({\bm{t}})\mathop{}\\!\mathrm{d}{\bm{t}}\right\\}}{\int_{\Omega}e^{-kf({\bm{t}})}\cdot\pi({\bm{t}})\mathop{}\\!\mathrm{d}{\bm{t}}}$
$\displaystyle\stackrel{{\scriptstyle}}{{=}}-f({\bm{x}})-\frac{\int_{\Omega}(-f({\bm{t}}))\cdot
e^{-kf({\bm{t}})}\cdot\pi({\bm{t}})d{\bm{t}}}{\int_{\Omega}e^{-kf({\bm{t}})}\cdot\pi({\bm{t}})\mathop{}\\!\mathrm{d}{\bm{t}}}$
$\displaystyle=\mathbb{E}_{m_{k}}(f)-f({\bm{x}}).$

Hence,

$\displaystyle\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}k}m_{k}({\bm{x}})=m_{k}({\bm{x}})\cdot\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}k}\log
m_{k}({\bm{x}})=m_{k}({\bm{x}})\cdot\left(\mathbb{E}_{m_{k}}(f)-f({\bm{x}})\right).$

Then we have

$\displaystyle m^{(k+\Delta k)}({\bm{x}})$
$\displaystyle=m_{k}({\bm{x}})+\int_{k}^{k+\Delta
k}\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}{\bm{v}}}m^{({\bm{v}})}({\bm{x}})\mathop{}\\!\mathrm{d}{\bm{v}}$
$\displaystyle=m_{k}({\bm{x}})+\int_{k}^{k+\Delta
k}m^{({\bm{v}})}({\bm{x}})\cdot\left[\mathbb{E}^{({\bm{v}})}(f)-f({\bm{x}})\right]\mathop{}\\!\mathrm{d}{\bm{v}}.$

Then there exist a $\xi\in(k,k+\Delta k)$ such that

$\displaystyle\frac{\mathbb{E}^{(k+\Delta k)}(f)-\mathbb{E}_{m_{k}}(f)}{\Delta
k}$ $\displaystyle=\frac{1}{\Delta
k}\int_{\Omega}f({\bm{x}})\left(m^{(k+\Delta
k)}({\bm{x}})-m_{k}({\bm{x}})\right)\mathrm{d}{\bm{x}}$
$\displaystyle=\frac{1}{\Delta k}\int_{\Omega}\int_{k}^{k+\Delta
k}f({\bm{x}})m^{({\bm{v}})}({\bm{x}})\left[\mathbb{E}^{({\bm{v}})}(f)-f({\bm{x}})\right]\mathrm{d}{\bm{v}}\mathrm{~{}d}{\bm{x}}$
$\displaystyle\stackrel{{\scriptstyle(R\ref{rmk:3})}}{{=}}\frac{1}{\Delta
k}\int_{k}^{k+\Delta
k}\int_{\Omega}f({\bm{x}})m^{({\bm{v}})}({\bm{x}})\left[\mathbb{E}^{({\bm{v}})}(f)-f({\bm{x}})\right]\mathrm{d}{\bm{x}}\mathrm{~{}d}{\bm{v}}$
$\displaystyle=\frac{1}{\Delta k}\int_{k}^{k+\Delta
k}\left\\{\left[\mathbb{E}^{({\bm{v}})}(f)\right]^{2}-\mathbb{E}^{({\bm{v}})}(f^{2})\right\\}\mathrm{d}{\bm{v}}$
$\displaystyle=\left[\mathbb{E}^{(\xi)}(f)\right]^{2}-\mathbb{E}^{(\xi)}(f^{2}),$

where the exchangeability of the integral is proved by Fubini Theorem in
Remark 3. Hence, we have

$\displaystyle\frac{\mathrm{d}\mathbb{E}_{m_{k}}(f)}{\mathrm{d}k}$
$\displaystyle=\lim_{\Delta k\rightarrow 0}\frac{\mathbb{E}^{(k+\Delta
k)}(f)-\mathbb{E}_{m_{k}}(f)}{\Delta k}$
$\displaystyle=\left\\{\mathbb{E}_{m_{k}}(f)\right\\}^{2}-\mathbb{E}_{(k)}(f^{2})$
$\displaystyle=-\int_{\Omega}\left(f({\bm{x}})-\mathbb{E}_{m_{k}}(f)\right)^{2}m_{k}({\bm{x}})\mathrm{d}{\bm{x}}$
$\displaystyle=-\mathbb{V}{ar}^{(k)}(f)\leqslant 0,$

where the equality holds only when $\mathbb{V}{ar}^{(k)}(f)=0$, i.e., $f$ is a
constant function on $\Omega$.

######  Remark 3 (R3).

Define $h_{k}(t)=te^{-kt}$ for any $t\in\mathbb{R}$, then its first derivative
is $h_{k}^{\prime}(t)=(1-kt)e^{-kt}$. For any $t<\frac{1}{k}$,
$h_{k}^{\prime}(t)>0$, and thus $g(t)$ is increasing when $t<\frac{1}{k}$; for
any $t\geqslant\frac{1}{k}$, $h_{k}^{\prime}(t)\leqslant 0$, and thus $g(t)$
is non-increasing when $t\geqslant\frac{1}{k}$. Therefore, $h_{k}(t)\leqslant
h_{k}(t=\frac{1}{k})=\frac{1}{ke}$. Define $g_{k}(t)=t^{2}e^{-kt}$, then the
absolute value of $g_{k}(t)$ has its upper bound as $\frac{4}{k^{2}e^{2}}$
using similar strategy. For any $k>0$, the denominator of
$\mathbb{E}^{({\bm{v}})}(f)$ is a finite positive constant, where
$0<\int_{\Omega}e^{-kf({\bm{x}})}\pi({\bm{x}})\mathrm{d}{\bm{x}}\stackrel{{\scriptstyle\Delta}}{{=}}\alpha_{k}\leqslant
e^{-kf^{*}}\int_{\Omega}\pi({\bm{x}})\mathrm{d}{\bm{x}}\leqslant e^{-kf^{*}}$,
and thus $\mathbb{E}^{({\bm{v}})}(f)$ defined in the proof of Theorem 2 is
bounded following:

$\displaystyle\mathbb{E}^{({\bm{v}})}(f)$
$\displaystyle=\frac{\int_{\Omega}f({\bm{x}})e^{-{\bm{v}}f({\bm{x}})}\pi({\bm{x}})\mathrm{d}{\bm{x}}}{\int_{\Omega}e^{-{\bm{v}}f({\bm{x}})}\pi({\bm{x}})\mathrm{d}{\bm{x}}}$
(12)
$\displaystyle=\frac{\int_{\Omega}h_{v}(f({\bm{x}}))\pi({\bm{x}})\mathrm{d}{\bm{x}}}{\int_{\Omega}e^{-{\bm{v}}f({\bm{x}})}\pi({\bm{x}})\mathrm{d}{\bm{x}}}$
$\displaystyle\leqslant\frac{\frac{1}{ve}\int_{\Omega}\pi({\bm{x}})\mathrm{d}{\bm{x}}}{\int_{\Omega}e^{-{\bm{v}}f({\bm{x}})}\pi({\bm{x}})\mathrm{d}{\bm{x}}}=\frac{1}{ve\alpha_{v}}$

Similarly, $\mathbb{E}^{({\bm{v}})}(f^{2})$ is also bounded by:

$\displaystyle\mathbb{E}^{({\bm{v}})}(f^{2})$
$\displaystyle=\frac{\int_{\Omega}f^{2}({\bm{x}})e^{-{\bm{v}}f({\bm{x}})}\pi({\bm{x}})\mathrm{d}{\bm{x}}}{\int_{\Omega}e^{-{\bm{v}}f({\bm{x}})}\pi({\bm{x}})\mathrm{d}{\bm{x}}}$
(13)
$\displaystyle=\frac{\int_{\Omega}g_{v}(f({\bm{x}}))\pi({\bm{x}})\mathrm{d}{\bm{x}}}{\int_{\Omega}e^{-{\bm{v}}f({\bm{x}})}\pi({\bm{x}})\mathrm{d}{\bm{x}}}$
$\displaystyle\leqslant\frac{\frac{4}{{\bm{v}}^{2}e^{2}}\int_{\Omega}\pi({\bm{x}})\mathrm{d}{\bm{x}}}{\int_{\Omega}e^{-{\bm{v}}f({\bm{x}})}\pi({\bm{x}})\mathrm{d}{\bm{x}}}=\frac{4}{{\bm{v}}^{2}e^{2}\alpha_{v}}$

Hence, given bounded $\mathbb{E}^{({\bm{v}})}(f)$ and bounded
$\mathbb{E}^{({\bm{v}})}(f)$, the following order of integration is
exchangeable by the Fubini theorem:

$\int_{\Omega}\int_{k}^{k+\Delta
k}\mathbb{E}^{({\bm{v}})}(f)\mathrm{d}{\bm{v}}\mathrm{~{}d}{\bm{x}}=\int_{k}^{k+\Delta
k}\int_{\Omega}\mathbb{E}^{({\bm{v}})}(f)\mathrm{d}{\bm{x}}\mathrm{~{}d}{\bm{v}},$
$\int_{\Omega}\int_{k}^{k+\Delta
k}\mathbb{E}^{({\bm{v}})}(f^{2})\mathrm{d}{\bm{v}}\mathrm{~{}d}{\bm{x}}=\int_{k}^{k+\Delta
k}\int_{\Omega}\mathbb{E}^{({\bm{v}})}(f^{2})\mathrm{d}{\bm{x}}\mathrm{~{}d}{\bm{v}}.$

∎

###  A.3 Proof for Theroem 3

###### Proof for Theroem 3.

First, we prove eq. 6.

For $\forall k>0$, denote
$c_{k}=\log\left\\{\int_{\Omega}e^{-kf({\bm{t}})}\pi({\bm{t}})\mathop{}\\!\mathrm{d}{\bm{t}}\right\\}$,
then

$\displaystyle\operatorname*{argmax}_{{\bm{x}}\in\Omega}m_{k}({\bm{x}})$
$\displaystyle=\operatorname*{argmax}_{{\bm{x}}\in\Omega}\log m_{k}({\bm{x}})$
$\displaystyle=\operatorname*{argmax}_{{\bm{x}}\in\Omega}\left(-kf({\bm{x}})+\log\pi({\bm{x}})-c_{k}\right)$
$\displaystyle=\operatorname*{argmin}_{{\bm{x}}\in\Omega}\left(f({\bm{x}})-\frac{1}{k}\log\pi({\bm{x}})\right),$

which indicates the maximizers of $m_{k}({\bm{x}})$ are the minimizers of
$\left(f({\bm{x}})-\frac{1}{k}\log\pi({\bm{x}})\right)$. Assume the upper
bound for $\pi({\bm{x}})$ is $B$, then for any ${\bm{x}}_{k}^{*}\in\Omega_{k}$
and $\tilde{{\bm{x}}}\in\Omega_{*}$, we have

$f({\bm{x}}_{k}^{*})-\frac{1}{k}\log B\leqslant
f({\bm{x}}_{k}^{*})-\frac{1}{k}\log\pi({\bm{x}}_{k}^{*})\leqslant
f(\tilde{{\bm{x}}})-\frac{1}{k}\log\pi(\tilde{{\bm{x}}}),$ (14)

where the first inequality of eq. (14) holds for the bounded density, where
$\pi({\bm{x}})\leqslant B,\forall{\bm{x}}\in\Omega$. The second inequality of
eq. (14) lies in ${\bm{x}}_{k}^{*}\in\arg\max m_{k}({\bm{x}})$. Furtherly,
$\lim_{k\rightarrow\infty}f({\bm{x}}_{k}^{*})=f^{*}$ follows from

$f^{*}\leqslant\liminf_{k\rightarrow\infty}f({\bm{x}}_{k}^{*})\leqslant\limsup_{k\rightarrow\infty}f({\bm{x}}_{k}^{*})\leqslant\limsup_{k\rightarrow\infty}\left\\{f(\tilde{{\bm{x}}})-\frac{1}{k}\log\pi(\tilde{{\bm{x}}})\right\\}=f^{*},$
(15)

where the three inequalities follow from the definition of $f^{*}$ ($\leqslant
f({\bm{x}})$, for $\forall{\bm{x}}\in\Omega$), the definition of limit
inferior and limit superior, and the limit superior of eq. (14) respectively.

Next, we prove eq. (7). In the previous part, we have already proved
$\lim_{k\rightarrow\infty}f({\bm{x}}_{k}^{*})=f^{*},$ i.e., for any
$\epsilon>0$, there exists large $K_{0}>0$, such that for $\forall k>K_{0}$,
it holds that

$f({\bm{x}}_{k}^{*})-f^{*}<\epsilon.$ (16)

Suppose there exists a large $K>K_{0}$, such that

$\inf_{\tilde{{\bm{x}}}\in\Omega^{*}}||{\bm{x}}_{K}^{*}-\tilde{{\bm{x}}}||>\delta,\quad\text{for
}\forall\tilde{{\bm{x}}}\in\Omega^{*}.$

Then $f({\bm{x}}_{K}^{*})-f^{*}>\epsilon$ given the strong separability
condition in Assumption 2, which is contradictory to eq. (16). Thus, we have
completed the proof for Theorem 3 using the proof by contradictory technique.
∎

######  Theorem 4.

The nascent minima distribution function defined in Definition 1 satisfies:

  1. (i)

For $\forall k\in\mathbb{R}$, $m_{k}({\bm{x}})$ is a PDF on $\Omega$;
especially when $k=0$, $m^{(0)}({\bm{x}})=\pi({\bm{x}})/\Pi(\Omega)$, where
$\Pi(S)=\int_{S}\pi({\bm{x}})\mathrm{d}{\bm{x}}$ is the probability measure of
$S\subseteq\mathbb{R}^{d}$.

  2. (ii)

If $\nabla\pi$ and $\nabla f:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ are
continuous real functions on each dimension, then

$\nabla\log m_{k}(x)=-k\nabla f(x)+\frac{\nabla\pi(x)}{\pi(x)}.$

  3. (ii)

For every $k\in\mathbb{R}$, it holds that

$\displaystyle\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}k}\log
m_{k}({\bm{x}})=\mathbb{E}_{m_{k}}(f)-f({\bm{x}}),$

where
$\mathbb{E}_{m_{k}}(f)=\int_{\Omega}f({\bm{x}})m_{k}({\bm{x}})\mathrm{d}{\bm{x}}$.

  4. (iii)

For any ${\bm{x}}\in\Omega$, consider the sequence $m_{k}({\bm{x}})$, where
$k>0$ and is going to infinity:

     * (a)

If $\Omega^{*}$ has zero probability measure, i.e.,
$\Pi(\Omega^{*})\triangleq\int_{\Omega^{*}}\pi({\bm{x}})\mathrm{d}{\bm{x}}=0$,
then

$\displaystyle
m_{\infty}({\bm{x}})=\lim_{k\rightarrow\infty}m_{k}({\bm{x}})=\left\\{\begin{array}[]{cc}0,&{\bm{x}}\notin\Omega^{*};\\\
\infty,&{\bm{x}}\in\Omega^{*}.\end{array}\right.$

     * (b)

If $\Omega^{*}$ has nonzero probability measure, i.e., $\Pi(\Omega^{*})>0$,
then

$\displaystyle\lim_{k\rightarrow\infty}m_{k}({\bm{x}})=\left\\{\begin{array}[]{cc}0,&{\bm{x}}\notin\Omega^{*};\\\
\frac{\pi({\bm{x}})}{\Pi(\Omega^{*})},&{\bm{x}}\in\Omega^{*}.\end{array}\right.$

  5. (iv)

Define $\Omega_{k}$ as the set of maximizers of $m_{k}({\bm{x}})$ on
${\bm{x}}\in\Omega$, if $\pi({\bm{x}})$ is bounded, then for any $k>0$, the
sequence of ${\bm{x}}_{k}^{*}\in\Omega_{k}$ satisfies the following:

    1. (a)

$\lim_{k\rightarrow\infty}f({\bm{x}}_{k}^{*})=0(=f^{*})$.

    2. (b)

$\pi({\bm{x}}_{k}^{*})\geqslant\pi^{*}\stackrel{{\scriptstyle\Delta}}{{=}}\max_{{\bm{x}}^{*}\in\Omega^{*}}\pi({\bm{x}}^{*})$.

    3. (c)

The sequence $\\{\pi({\bm{x}}_{k}^{*})\\}_{k=1,2,\cdots}$ is non-increasing
and converges to a limit, where
$\lim_{k\rightarrow\infty}\pi({\bm{x}}_{k}^{*})=\liminf_{k\rightarrow\infty}\pi({\bm{x}}_{k}^{*})\stackrel{{\scriptstyle\Delta}}{{=}}\pi_{0}$.

###### Proof.

Clearly, (i) follows from the definition of $m_{k}({\bm{x}})$ in Definition 1.

  * •

For every $k\in\mathbb{R}$, (ii) follows from

$\displaystyle\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}k}\log
m_{k}({\bm{x}})$
$\displaystyle=\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}k}\left(-kf({\bm{x}})+\log\pi({\bm{x}})-\log\left\\{\int_{\Omega}e^{-kf(t)}\cdot\pi(t)dt\right\\}\right)$
$\displaystyle\stackrel{{\scriptstyle}}{{=}}-f({\bm{x}})-\frac{\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}k}\left\\{\int_{\Omega}e^{-kf(t)}\cdot\pi(t)dt\right\\}}{\int_{\Omega}e^{-kf(s)}\cdot\pi(s)ds}$
$\displaystyle\stackrel{{\scriptstyle}}{{=}}-f({\bm{x}})-\frac{\int_{\Omega}(-f(t))\cdot
e^{-kf(t)}\cdot\pi(t)dt}{\int_{\Omega}e^{-kf(s)}\cdot\pi(s)ds}$
$\displaystyle=\mathbb{E}_{m_{k}}(f)-f({\bm{x}}),$

where
$\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}k}\left\\{\int_{\Omega}e^{-kf(t)}\cdot\pi(t)dt\right\\}=\int_{\Omega}\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}k}\left\\{e^{-kf(t)}\cdot\pi(t)\right\\}dt$
follows from the Tonelli theorem for the exchangeable order of the integral
and derivatives, given that the function
$\left\\{e^{-kf(t)}\cdot\pi(t)\right\\}$ is non-negative.

  * •

For (iii) notice that by Remark 1,

$\int_{\Omega}e^{-kf({\bm{x}})}\pi({\bm{x}})\mathrm{d}{\bm{x}}=\int_{\Omega^{*}}\pi({\bm{x}})\mathrm{d}{\bm{x}}+\int_{\Omega\setminus\Omega^{*}}e^{-kf({\bm{x}})}\pi({\bm{x}})\mathrm{d}{\bm{x}},$

since $f({\bm{x}})=0$ for $\forall{\bm{x}}\in\Omega^{*}$. In addition, since
$\left\\{e^{-kf({\bm{x}})}\pi({\bm{x}})\right\\}$ is monotonely decreasing as
$k$ increases and $\lim_{k\rightarrow\infty
0}\left\\{e^{-kf({\bm{x}})}\pi({\bm{x}})\right\\}=0$ due to $f({\bm{x}})>0$
for any ${\bm{x}}\notin\Omega^{*}$, then

$\lim_{k\rightarrow\infty}\int_{\Omega}e^{-kf({\bm{x}})}\pi({\bm{x}})\mathrm{d}{\bm{x}}=\int_{\Omega^{*}}\pi({\bm{x}})\mathrm{d}{\bm{x}}+\lim_{k\rightarrow\infty}\int_{\Omega\setminus\Omega^{*}}e^{-kf({\bm{x}})}\pi({\bm{x}})\mathrm{d}{\bm{x}}=\Pi(\Omega^{*}),$

which follows from the monotone convergence theorem.

    * (a)

Hence, for any ${\bm{x}}\in\Omega^{*}$,

$m_{\infty}({\bm{x}})=\lim_{k\rightarrow\infty}m_{k}({\bm{x}})=\lim_{k\rightarrow\infty}\frac{e^{-kf({\bm{x}})}\cdot\pi({\bm{x}})}{\int_{\Omega}e^{-kf(t)}\cdot\pi(t)dt}=\left\\{\begin{array}[]{cc}\infty,&\text{if
}\Pi(\Omega^{*})=0;\\\ \frac{\pi({\bm{x}})}{\Pi(\Omega^{*})},&\text{if
}\Pi(\Omega^{*})\neq 0.\end{array}\right.$

    * (b)

For any ${\bm{x}}^{\prime}\notin\Omega^{*}$, it follows from the continuity of
$f$ that there exists a set $\Omega_{{\bm{x}}^{\prime}}$ such that
$f(t)<f({\bm{x}}^{\prime})$ for any $t\in\Omega_{{\bm{x}}^{\prime}}$, hence,

$\displaystyle m_{k}({\bm{x}}^{\prime})$
$\displaystyle=\frac{e^{-kf({\bm{x}}^{\prime})}\cdot\pi({\bm{x}}^{\prime})}{\int_{\Omega_{{\bm{x}}^{\prime}}}e^{-kf(t)}\cdot\pi(t)dt+\int_{\Omega\setminus\Omega_{{\bm{x}}^{\prime}}}e^{-kf(t)}\cdot\pi(t)dt}$
$\displaystyle\leqslant\frac{\pi({\bm{x}}^{\prime})}{\int_{\Omega_{{\bm{x}}^{\prime}}}(e^{-f(t)}/e^{-f({\bm{x}}^{\prime})})^{k}\cdot\pi(t)dt},$

since $e^{-f(t)}/e^{-f({\bm{x}}^{\prime})}>1$ for any
$t\in\Omega_{{\bm{x}}^{\prime}}$ and $\pi(t)>0$, the limit of
$\int_{\Omega_{{\bm{x}}^{\prime}}}(e^{-f(t)}/e^{-f({\bm{x}}^{\prime})})^{k}\cdot\pi(t)dt$
tends to $\infty$ as $k\rightarrow\infty$; thus, it holds that

$\displaystyle\operatorname*{lim}_{k\rightarrow\infty}m_{k}({\bm{x}}^{\prime})=0,\forall{\bm{x}}^{\prime}\notin\Omega^{*}.$

  * •

(iv) describes the properties of the maximizers of $m_{k}({\bm{x}})$. For
$\forall k>0$, denote
$c_{k}=\log\left\\{\int_{\Omega}e^{-kf(t)}\pi(t)dt\right\\}$, then

$\displaystyle\operatorname*{argmax}_{{\bm{x}}\in\Omega}m_{k}({\bm{x}})$
$\displaystyle=\operatorname*{argmax}_{{\bm{x}}\in\Omega}\log m_{k}({\bm{x}})$
$\displaystyle=\operatorname*{argmax}_{{\bm{x}}\in\Omega}\left(-kf({\bm{x}})+\log\pi({\bm{x}})-c_{k}\right)$
$\displaystyle=\operatorname*{argmin}_{{\bm{x}}\in\Omega}\left(f({\bm{x}})-\frac{1}{k}\log\pi({\bm{x}})\right)$
$\displaystyle=\operatorname*{argmin}_{{\bm{x}}\in\Omega}g_{k}({\bm{x}}),$

which indicates the maximizers of $m_{k}({\bm{x}})$ are the minimizers of
$g_{k}({\bm{x}})$, i.e.,
$\Omega_{k}=\\{{\bm{x}}^{\prime}:m_{k}({\bm{x}}^{\prime})\geqslant
m_{k}({\bm{x}}),\forall{\bm{x}}\in\Omega\\}=\\{{\bm{x}}^{\prime\prime}:g_{k}({\bm{x}}^{\prime\prime})\leqslant
g_{k}({\bm{x}}),\forall{\bm{x}}\in\Omega\\}$.

    1. (a)

For any ${\bm{x}}_{k}^{*}\in\Omega_{k}$ and ${\bm{x}}^{*}\in\Omega^{*}$, we
have

$f({\bm{x}}_{k}^{*})-\frac{1}{k}\log B\leqslant
f({\bm{x}}_{k}^{*})-\frac{1}{k}\log\pi({\bm{x}}_{k}^{*})\leqslant
f({\bm{x}}^{*})-\frac{1}{k}\log\pi({\bm{x}}^{*}),$

where the first inequality lies in $\pi({\bm{x}})\leqslant
B,\forall{\bm{x}}\in\Omega$ and the second inequality lies in
${\bm{x}}_{k}^{*}\in\arg\max m_{k}({\bm{x}})$. Furtherly,
$\lim_{k\rightarrow\infty}f({\bm{x}}_{k}^{*})=f^{*}$ follows from

$f^{*}\leqslant\liminf_{k\rightarrow\infty}f({\bm{x}}_{k}^{*})\leqslant\limsup_{k\rightarrow\infty}f({\bm{x}}_{k}^{*})\leqslant\limsup_{k\rightarrow\infty}\left\\{f({\bm{x}}^{*})-\frac{1}{k}\log\pi({\bm{x}}^{*})\right\\}=f^{*}.$

    2. (b)

For any ${\bm{x}}^{*}\in\Omega^{*}$, and ${\bm{x}}_{k}^{*}\in\Omega_{k},$ we
have
$\pi({\bm{x}}_{k}^{*})\geqslant\pi^{*}=\max_{{\bm{x}}\in\Omega^{*}}\pi({\bm{x}}^{*})$
for $k>0$, which follows from:

$f({\bm{x}}_{k}^{*})-\frac{1}{k}\log\pi({\bm{x}}_{k}^{*})\leqslant
f({\bm{x}}^{*})-\frac{1}{k}\log\pi({\bm{x}}^{*})\leqslant
f({\bm{x}}_{k}^{*})-\frac{1}{k}\log\pi({\bm{x}}^{*}).$

    3. (c)

The monotonicity of the sequence $\\{\pi({\bm{x}}_{k}^{*})\\}_{k=1,2,\cdots}$
follows from:

$\displaystyle f({\bm{x}}_{k+1}^{*})-\frac{1}{k+1}\log\pi({\bm{x}}_{k+1}^{*})$
$\displaystyle\leqslant
f({\bm{x}}_{k}^{*})-\frac{1}{k+1}\log\pi({\bm{x}}_{k}^{*})$ $\displaystyle\iff
f({\bm{x}}_{k+1}^{*})-\frac{1}{k+1}\log\pi({\bm{x}}_{k+1}^{*})$
$\displaystyle\leqslant
f({\bm{x}}_{k}^{*})-\frac{1}{k}\log\pi({\bm{x}}_{k}^{*})+\frac{1}{k(k+1)}\log\pi({\bm{x}}_{k}^{*})$
$\displaystyle\iff
f({\bm{x}}_{k+1}^{*})-\frac{1}{k+1}\log\pi({\bm{x}}_{k+1}^{*})$
$\displaystyle\leqslant
f({\bm{x}}_{k+1}^{*})-\frac{1}{k}\log\pi({\bm{x}}_{k+1}^{*})+\frac{1}{k(k+1)}\log\pi({\bm{x}}_{k}^{*})$
$\displaystyle\iff\log\pi({\bm{x}}_{k+1}^{*})$
$\displaystyle\leqslant\log\pi({\bm{x}}_{k}^{*})$

∎

##  Appendix B Experimental Settings

Consistent with the experimental parameters adopted by ZoRD (Shu et al. 2022),
we designate the domains of the Ackley and Levy functions as $[-20,20]^{d}$
and $[-7.5,7.5]^{d}$, respectively. For ProGO, the parameter for LSS is sample
size as $N=200$ and burn-in period $n_{b}=20$. The configurations are
uniformly applied across the RGF, PRGF, GD, and ZoRD algorithms to ensure fair
comparisons, where the Adam optimizer (Kingma & Ba 2015) is utilized with a
fixed learning rate of 0.1 and exponential decay rates of 0.9 and 0.999. It
should be noted that while ProGO is implemented in the R environment, other
algorithms are executed in Python. Although this discrepancy may affect
runtime comparisons, it does not influence accuracy results.