# Simulated annealing from continuum to discretization: a convergence analysis
via the Eyring–Kramers law

Wenpin Tang  Department of Industrial Engineering and Operations Research,
Columbia University<EMAIL_ADDRESS>and  Xun Yu Zhou  Department of
Industrial Engineering and Operations Research, Columbia University.
<EMAIL_ADDRESS>
###### Abstract.

We study the convergence rate of continuous-time simulated annealing
$(X_{t};\,t\geq 0)$ and its discretization $(x_{k};\,k=0,1,\ldots)$ for
approximating the global optimum of a given function $f$. We prove that the
tail probability $\mathbb{P}(f(X_{t})>\min f+\delta)$ (resp.
$\mathbb{P}(f(x_{k})>\min f+\delta)$) decays polynomial in time (resp. in
cumulative step size), and provide an explicit rate as a function of the model
parameters. Our argument applies the recent development on functional
inequalities for the Gibbs measure at low temperatures – the Eyring-Kramers
law. In the discrete setting, we obtain a condition on the step size to ensure
the convergence.

Key words: Simulated annealing, convergence rate, Euler discretization,
Eyring-Kramers law, functional inequalities, overdamped Langevin equation.

##  1\. Introduction

Simulated annealing (SA) is an umbrella term which denotes a set of stochastic
optimization methods. The goal of SA is to find the global minimum of a
function $f:\mathbb{R}^{d}\to\mathbb{R}$, in particular when $f$ is non-
convex. These methods have many applications in physics, operations research
and computer science, see e.g. van Laarhoven and Aarts (1987); Koulamas et al.
(1994); Delahaye et al. (2019). The name is inspired from annealing in
metallurgy, which is a process aiming to increase the size of crystals by
heating and controlled cooling. The stochastic version of SA was independently
proposed by Kirkpatrick et al. (1983) and Cerny (1985). The idea is as
follows: consider a stochastic process related to $f$ which is subject to
thermal noise. When simulating this process, one decreases the temperature
slowly over the time. As a result, the stochastic process escapes from saddle
points and local optima, and converges to the global minimum of $f$ with high
probability. This is generally true if the cooling is slow enough, and it is
important to find the right stochastic process with the fastest possible
cooling schedule which approximates the global optimum.

In this paper, we explore the convergence rate of continuous-time SA and its
discretization. To be more precise, define the continuous-time SA process
$(X_{t};\,t\geq 0)$ by

$dX_{t}=-\nabla f(X_{t})dt+\sqrt{2\tau_{t}}\,dB_{t},\quad
X_{0}\stackrel{{\scriptstyle d}}{{=}}\mu_{0}(dx),$ (1)

where $(B_{t};\,t\geq 0)$ is $d$-dimensional Brownian motion, $\tau_{t}$ is
the cooling schedule of temperature, and $\mu_{0}(dx)$ is some initial
distribution. This formulation was first considered by Geman and Hwang (1986).
If $\tau_{t}\equiv\tau$ is constant in time, the process (1) is the well-known
overdamped Langevin equation whose stationary distribution is the Gibbs
measure $\nu_{\tau}(dx)\propto\exp(-f(x)/\tau)dx$. Thus, we sometimes call the
process (1) an SA adapted overdamped Langevin equation as well. See Section 3
for background. We also consider the Euler-Maruyama discretization of the
continuous-time simulated annealing process. Let $\eta_{k}$ be the step size
at iteration $k$, and $\Theta_{k}:=\sum_{j\leq k}\eta_{j}$ be the cumulative
step size up to iteration $k$. The discrete SA process
$(x_{k};\,k=0,1,\ldots)$ is defined by

$x_{k+1}=x_{k}-\nabla
f(x_{k})\eta_{k}+\sqrt{2\tau_{\scaleto{\Theta_{k}}{5pt}}\eta_{k}}Z_{k},\quad
x_{0}\stackrel{{\scriptstyle d}}{{=}}\mu_{0}(dx),$ (2)

where $(Z_{k};\,k=0,1,\ldots)$ are independent and identically distributed
standard normal vectors, and $\tau_{\scaleto{\Theta_{k}}{5pt}}$ is the cooling
schedule at iteration $k$. For $\tau_{t}\equiv\tau$ constant in time, the
scheme (2) is known as the unadjusted Langevin algorithm (ULA) which
approximates the Gibbs measure $\nu_{\tau}$. The algorithm was introduced in
Parisi (1981); Grenander and Miller (1994), and further studied by Roberts and
Tweedie (1996); Durmus and Moulines (2017).

The goal of this paper is to study the decay in time of the tail probability

$\mathbb{P}(f(X_{t})>\min f+\delta)\quad\mbox{or}\quad\mathbb{P}(f(x_{k})>\min
f+\delta),$

under suitable conditions on the function $f$, the cooling schedule
$\tau_{t}$, and the discretization scheme $\eta_{k},\Theta_{k}$. Let us
mention a few motivations. First there are a line of works on the interplay
between sampling and optimization (Raginsky et al. (2017); Ma et al. (2019a,
b)). Note that if $\tau_{t}\equiv\tau$ is constant in time, the overdamped
Langevin equation converges to the Gibbs measure $\nu_{\tau}$, and for $\tau$
sufficiently small, the Gibbs measure $\nu_{\tau}$ approximates the Dirac mass
at the global minimum of $f$. Combining these two ingredients, SA sets the
cooling schedule $\tau_{t}$ to decrease to $0$ over the time. It is then
expected that the SA process (1) or (2) converges to the global minimum as
$t\to\infty$ or $k\to\infty$. There are also recent works on various noisy
gradient-based algorithms (Ge et al. (2015); Jin et al. (2017); Chen et al.
(2020); Guo et al. (2020)), which aims to escape saddle points and find a
local minimum of $f$. In comparison with these methods, SA has priority to
find the global minimum but at the cost of longer exploration time.

The main tool in our analysis is the Eyring-Kramers law, which is a set of
functional inequalities for the Gibbs measure at low temperatures. Let us
explain our results. It was shown in Geman and Hwang (1986); Chiang et al.
(1987) that the correct order of $\tau_{t}$ for the process (1) to converge to
the global minimum of $f$ is $(\ln t)^{-1}$. In fact, there is a phase
transition related to the critical depth $E_{*}$ of the function $f$:

  1. (a)

If $\tau_{t}\leq\frac{E}{\ln t}$ for $t$ large enough with $E<E_{*}$, then
$\limsup_{t\to\infty}\mathbb{P}(f(X_{t})\leq\min f+\delta)<1$.

  2. (b)

If $\tau_{t}\geq\frac{E}{\ln t}$ for $t$ large enough with $E>E_{*}$, then
$\lim_{t\to\infty}\mathbb{P}(f(X_{t})\leq\min f+\delta)=1$.

Roughly speaking, the critical depth $E_{*}$ is the largest hill one needs to
climb starting from a local minimum to the global minimum. The formal
definition of the critical depth $E_{*}$ will be given in Assumption 2, but
see Figure 1 below for an illustration when $f$ is a double-well function.

Figure 1. Illustration of the critical depth of a double-well function.

Part $(a)$ was proved by Holley et al. (1989) using a sophisticated argument
that involves the Poincaré inequality. Part $(b)$ was proved by Miclo (1992)
who characterized the fastest cooling schedule by the Eying-Kramers law for
the log-Sobolev inequality. See also Miclo (1995); Zitt (2008); Fournier and
Tardif (2020) for similar results under different conditions on the function
$f$. The convergence rate of $(f(X_{t});\,t\geq 0)$ to the global minimum of
$f$ was considered by Márquez (1997) only for $\delta$ sufficiently small.
However, the rate of $\mathbb{P}(f(X_{t})>\min f+\delta)$ for general
$\delta>0$ has yet been studied since no estimates of the log-Sobolev
inequality for the Gibbs measure at low temperatures were known until the
mid-2010s. Taking advantage of recently developed theory (Menz and Schlichting
(2014); Menz et al. (2018); Vempala and Wibisono (2019)), we are able to give
a non-asymptotic convergence rate of both continuous-time and discrete SA.

To simplify the notation, we assume throughout this paper that

$\min_{\mathbb{R}^{d}}f(x)=0,$

i.e. the global minimum of $f$ is $0$ by considering $f-\min f$. Our results
are outlined as follows, and the precise statement of these results will be
given in Section 2.

######  Main Results (Informal).

Under some assumptions on the function $f$, we have:

  1. (i)

Assume that $\tau_{t}\geq\frac{E}{\ln t}$ for $t$ large enough with $E>E_{*}$.
For $\delta>0$, there exists $C>0$ (depending on $\delta,f,E,d$) such that

$\mathbb{P}(f(X_{t})>\delta)\leq
Ct^{-\min\left(\frac{\delta}{E},\frac{1}{2}(1-\frac{E_{*}}{E})\right)}.$

  2. (ii)

Assume that $\tau_{t}\geq\frac{E}{\ln t}$ for $t$ large enough with $E>E_{*}$.
Also assume that $\Theta_{k}\to\infty$, $\eta_{k+1}\Theta_{k}\to 0$ and
$\frac{\ln\Theta_{k}}{\sum_{j\leq k}\eta_{j+1}\Theta_{j}^{-E_{*}/E}}\to 0$ as
$k\to\infty$. For $\delta>0$, there exists $C>0$ (depending on $\delta,f,E,d$)
such that

$\mathbb{P}(f(x_{k})>\delta)\leq
C\Theta_{k}^{-\min\left(\frac{\delta}{E},\frac{1}{2}(1-\frac{E_{*}}{E})\right)}.$

The contribution of this paper is twofold:

  * •

Polynomial decay rate. We prove that the tail probability
$\mathbb{P}(f(X_{t})>\delta)$ (resp. $\mathbb{P}(f(x_{k})>\delta)$) decays
polynomially in time (resp. in cumulative step size). In the continuous
setting, Monmarché (2018) obtained the same rate of convergence for SA adapted
underdamped Langevin equation, and Menz et al. (2018) considered an
improvement of SA via parallel tempering. However, the convergence rate for
continuous-time SA, i.e. part $(i)$ has yet been recorded in the literature
though this result is probably understood in the Eyring-Kramers folklore. We
provide a self-contained treatment which bridges the literature and, more
importantly, can be applied to obtain the new result in the discrete setting.

  * •

Choice of step size. Part $(ii)$ for the discrete simulated annealing is
completely novel to our best knowledge, and it also gives a practical guidance
on the choice of step size for discretization. The condition
$\Theta_{k}\to\infty$ indicates that the step size cannot be chosen too small,
while the condition $\eta_{k+1}\Theta_{k}\to 0$ suggests that the step size
cannot be chosen too large. The condition $\frac{\ln\Theta_{k}}{\sum_{j\leq
k}\eta_{j+1}\Theta_{j}^{-E_{*}/E}}\to 0$ is purely technical. For instance,
$\eta_{k}=k^{-\theta}$ with $\theta\in(\frac{1}{2},1]$ satisfies the
conditions in $(ii)$ to ensure the convergence.

Also note that the rate
$\min\left(\frac{\delta}{E},\frac{1}{2}(1-\frac{E_{*}}{E})\right)$ is smaller
than $\frac{1}{2}$. It is interesting to know whether this rate is optimal,
and we leave the problem for future work.

There is another interesting problem: the dependence of the constant $C$ on
the dimension $d$. The issue is subtle, since most analysis including the
Eyring-Kramers law uses Laplace’s method. However, Laplace’s method may fail
if both the dimension $d$ and the inverse temperature $1/\tau$ tend to
infinity (Shun and McCullagh (1995)). As shown in Remark 1, we obtain an upper
bound for $C$ which is exponential in $d$. This suggests the convergence rate
is exponentially slow as the dimension increases, which concurs with the fact
that finding the global minimum of a general nonconvex function is NP-hard
(Jain and Kar (2017)).

Finally, we mention a few approaches to accelerate or improve SA. Fang et al.
(1997) considered a cooling schedule depending on both time and state;
Pavlyukevich (2007) used the Lévy flight; Monmarché (2018) studied SA adapted
to underdamped Langevin equation; Menz et al. (2018) applied the replica
exchange technique; Gao et al. (2020) employed a relaxed stochastic control
formulation, originally proposed by Wang et al. (2020) for reinforcement
learning, to derive a state-dependent temperature control schedule.

The remainder of the paper is organized as follows. Section 2 presents the
assumptions and our main results. Section 3 provides background on diffusion
processes and functional inequalities. The result for the continuous-time
simulated annealing (Theorem 1) is proved in Section 4. The result for the
discrete simulated annealing (Theorem 2) is proved in Section 5. We conclude
with Section 6.

##  2\. Main results

In this section, we make precise the informal statement in the introduction,
and present the main results of the paper. We first collect the notations that
will be used throughout this paper.

  * –

The notation $|\cdot|$ is the Euclidean norm of a vector, and $a\cdot b$ is
the scalar product of vectors $a$ and $b$.

  * –

For a function $f:\mathbb{R}^{d}\to\mathbb{R}$, let $\nabla f$, $\nabla^{2}f$
and $\Delta f$ denote its gradient, Hessian and Laplacian respectively.

  * –

For $X,Y$ two random variables, $||X-Y||_{TV}$ denotes the total variation
norm of the signed measure corresponding to $X-Y$.

  * –

The symbol $a\sim b$ means that $a/b\to 1$ as some problem parameter tends to
$0$ or $\infty$. Similarly, the symbol $a=\mathcal{O}(b)$ means that $a/b$ is
bounded as some problem parameter tends to $0$ or $\infty$.

We use $C$ for a generic constant which depends on problem parameters
($\delta,f,E\ldots$), and may change from line to line.

Next, we present a few assumptions on the function $f$. These assumptions are
standard in the study of metastability.

######  Assumption 1.

Let $f:\mathbb{R}^{d}\to\mathbb{R}$ be smooth, bounded from below, and satisfy
the conditions:

  1. (i)

$f$ is non-degenerate on the set of critical points. That is, for some $C>0$,

$\frac{|\xi|}{C}\leq|\nabla^{2}f(x)\xi|\leq C|\xi|\quad\mbox{for each
}x\in\\{z:\nabla f(z)=0\\}.$

  2. (ii)

There exists $C>0$ such that

$\liminf_{|x|\to\infty}\frac{|\nabla f(x)|^{2}-\Delta f(x)}{|x|^{2}}\geq
C,\quad\inf_{x}\nabla^{2}f(x)\geq-C.$

Let us make a few comments on Assumption 1. The condition $(ii)$ implies that
$f$ has at least quadratic growth at infinity. This is a necessary and
sufficient condition to obtain the log-Sobolev inequality (see (Royer, 2007,
Theorem 3.1.21) and Section 3.2) which is key to convergence analysis. The
conditions $(i)$ and $(ii)$ imply that the set of critical points is discrete
and finite (Menz and Schlichting, 2014, Remark 1.6). In particular, it follows
that the set of local minimum points $\\{m_{1},\ldots,m_{N}\\}$ is also
finite, with $N$ the number of local minimum points of $f$.

To keep the presentation simple, we make additional assumptions on $f$,
following (Menz et al., 2018, Assumption 2.5). Define the saddle height
$\widehat{f}(m_{i},m_{j})$ between two local minimum points $m_{i},m_{j}$ by

$\widehat{f}(m_{i},m_{j}):=\inf\left\\{\max_{s\in[0,1]}f(\gamma(s)):\gamma\in\mathcal{C}[0,1],\,\gamma(0)=m_{i},\,\gamma(1)=m_{j}\right\\}.$
(3)

See Figure 1 for an illustration of the saddle height
$\widehat{f}(m_{0},m_{1})$ when $f$ is a double-well function with $m_{0}$ the
global minimum and $m_{1}$ the local minimum.

######  Assumption 2.

Let $m_{1},\cdots,m_{N}$ be the positions of the local minima of $f$.

  1. (i)

$m_{1}$ is the unique global minimum point of $f$, and $m_{1},\ldots,m_{N}$
are ordered in the sense that there exists $\delta>0$ such that

$f(m_{N})\geq f(m_{N-1})\geq\cdots\geq f(m_{2})\geq\delta\quad\mbox{and}\quad
f(m_{1})=0.$

  2. (ii)

For each $i,j\in\\{1,\ldots,N\\}$, the saddle height between $m_{i},m_{j}$ is
attained at a unique critical point $s_{ij}$ of index one. That is,
$f(s_{ij})=\widehat{f}(m_{i},m_{j})$, and if
$\\{\lambda_{1},\ldots,\lambda_{n}\\}$ are the eigenvalues of
$\nabla^{2}f(s_{ij})$, then $\lambda_{1}<0$ and $\lambda_{i}>0$ for
$i\in\\{2,\ldots,n\\}$. The point $s_{ij}$ is called the communicating saddle
point between the minima $m_{i}$ and $m_{j}$.

  3. (iii)

There exists $p\in[N]$ such that the energy barrier $f(s_{p1})-f(m_{p})$
dominates all the others. That is, there exists $\delta>0$ such that for all
$i\in[N]\setminus\\{p\\}$,

$E_{*}:=f(s_{p1})-f(m_{p})\geq f(s_{i1})-f(m_{i})+\delta.$

The dominating energy barrier $E_{*}$ is called the critical depth.

The convergence result for the continuous-time SA (1) is stated as follows.
The proof will be given in Section 4.

######  Theorem 1.

Let $f$ satisfy Assumption 1 $\&$ 2. Assume that $\tau_{t}\sim\frac{E}{\ln t}$
and
$\frac{d}{dt}\left(\frac{1}{\tau_{t}}\right)=\mathcal{O}\left(\frac{1}{t}\right)$
as $t\to\infty$ with $E>E_{*}$. Also assume the moment condition for the
initial distribution $\mu_{0}$: for each $p\geq 1$, there exists $C_{p}>0$
such that

$\int_{\mathbb{R}^{d}}f(x)^{p}\mu_{0}(dx)\leq C_{p}.$ (4)

Then for each $\delta,\varepsilon>0$, there exists $C>0$ (depending on
$\delta,\varepsilon,f,\mu_{0},E,d$) such that

$\mathbb{P}(f(X_{t})>\delta)\leq
Ct^{-\min\left(\frac{\delta}{E},\frac{1}{2}(1-\frac{E_{*}}{E})\right)+\varepsilon}.$
(5)

To get the convergence result for the discrete simulated annealing, we need an
additional condition on the function $f$.

######  Assumption 3.

The gradient $\nabla f$ is globally Lipschitz, i.e. $|\nabla f(x)-\nabla
f(y)|\leq L|x-y|$ for some $L>0$.

The convergence result for the discrete simulated annealing (2) is stated as
follows. The proof will be given in Section 5.

######  Theorem 2.

Let $f$ satisfy Assumption 1, 2 $\&$ 3, and let the condition (4) for
$\mu_{0}$ hold. Assume that $\tau_{t}\sim\frac{E}{\ln t}$ and
$\frac{d}{dt}\left(\frac{1}{\tau_{t}}\right)=\mathcal{O}\left(\frac{1}{t}\right)$
as $t\to\infty$ with $E>E_{*}$. Also assume that $\Theta_{k}\to\infty$,

$\eta_{k+1}\Theta_{k}\to 0,$ (6)

and

$\frac{\ln\Theta_{k}}{\sum_{j\leq k}\eta_{j+1}\Theta_{j}^{-E_{*}/E}}\to 0,$
(7)

as $k\to\infty$. Then for each $\delta,\varepsilon>0$, there exists $C>0$
(depending on $\delta,\varepsilon,f,\mu_{0},E,d$) such that

$\mathbb{P}(f(x_{k})>\delta)\leq
C\Theta_{k}^{-\min\left(\frac{\delta}{E},\frac{1}{2}(1-\frac{E_{*}}{E})\right)+\varepsilon}.$
(8)

##  3\. Preliminaries

In this section, we present a few vocabularies and basic results of diffusion
processes and functional inequalities. We also explain how these results are
applied in the setting of SA, which will be useful in our convergence
analysis.

###  3.1. Diffusion processes and SA

Consider the general diffusion process $(X_{t};\,t\geq 0)$ in $\mathbb{R}^{d}$
of form:

$dX_{t}=b(t,X_{t})dt+\sigma(t,X_{t})dB_{t},\quad X_{0}\stackrel{{\scriptstyle
d}}{{=}}\mu_{0}(dx),$ (9)

where $(B_{t};t\geq 0)$ is a $d$-dimensional Brownian motion, with the drift
$b:\mathbb{R}_{+}\times\mathbb{R}^{d}\to\mathbb{R}^{d}$ and the covariance
matrix $\sigma:\mathbb{R}_{+}\times\mathbb{R}^{d}\to\mathbb{R}^{d\times d}$.
To ensure the well-posedness of the SDE (9), it requires some growth and
regularity conditions on $b$ and $\sigma$. For instance,

  * •

If $b$ and $\sigma$ are Lipschitz and have linear growth in $x$ uniformly in
$t$, then (9) has a strong solution which is pathwise unique.

  * •

If $b$ is bounded, and $\sigma$ is bounded, continuous and strictly elliptic,
then (9) has a weak solution which is unique in distribution.

We refer to Stroock and Varadhan (1979); Rogers and Williams (1987) for
background and further developments on the well-posedness of SDEs, and to
(Cherny and Engelbert, 2005, Chapter 1) for a review of related results.

Another important aspect is the distributional property of $(X_{t};\,t\geq 0)$
governed by the SDE (9). Let $\mathcal{L}$ be the infinitesimal generator of
the diffusion process $X$ defined by

$\displaystyle\mathcal{L}g(t,x)$ $\displaystyle:=b(t,x)\cdot\nabla
g(x)+\frac{1}{2}\sigma(t,x)\sigma(t,x)^{T}:\nabla^{2}g(x)$
$\displaystyle:=\sum_{i=1}^{d}b_{i}(t,x)\frac{\partial}{\partial
x_{i}}g(x)+\frac{1}{2}\sum_{i,j=1}^{d}\left(\sigma(t,x)\sigma(t,x)^{T}\right)_{ij}\frac{\partial^{2}}{\partial
x_{i}\partial x_{j}}g(x),$ (10)

and $\mathcal{L}^{*}$ be the corresponding adjoint operator given by

$\displaystyle\mathcal{L}^{*}g(t,x)$
$\displaystyle:=-\nabla\cdot(b(t,x)g(x))+\frac{1}{2}\nabla^{2}:(\sigma(t,x)\sigma(t,x)^{T}g(x))$
$\displaystyle:=-\sum_{i=1}^{d}\frac{\partial}{\partial
x_{i}}(b_{i}(t,x)g(x))+\frac{1}{2}\sum_{i,j=1}^{d}\frac{\partial^{2}}{\partial
x_{i}\partial x_{j}}(\sigma(t,x)\sigma(t,x)^{T}g(x))_{ij},$ (11)

where $g:\mathbb{R}^{d}\to\mathbb{R}$ is a suitably smooth test function, and
$:$ denotes the Frobenius inner product which is the component-wise inner
product of two matrices. The probability density $\rho_{t}(\cdot)$ of the
process $X$ at time $t$ then satisfies the Fokker-Planck equation:

$\frac{\partial\rho_{t}}{\partial t}=\mathcal{L}^{*}\rho_{t}.$ (12)

Specializing (12) to the SA process (1) with $b(t,x)=-\nabla f(x)$ and
$\sigma(t,x)=\sqrt{2\tau_{t}}\,I_{d}$, we have that the probability density
$\mu_{t}(\cdot)$ of $X$ governed by the SDE (1) satisfies

$\frac{\partial\mu_{t}}{\partial t}=\nabla\cdot(\mu_{t}\nabla
f)+\tau_{t}\Delta\mu_{t}.$ (13)

Under further growth conditions on $b$ and $\sigma$, it can be shown that as
$t\to\infty$, $\rho_{t}(\cdot)\to\rho_{\infty}(\cdot)$ which is the stationary
distribution of $(X_{t};\,t\geq 0)$. It is easily deduced from (12) that
$\rho_{\infty}$ is characterized by the equation
$\mathcal{L}^{*}\rho_{\infty}=0$; see Ethier and Kurtz (1986); Meyn and
Tweedie (1993) for a general theory on stability of diffusion processes, and
(Tang, 2019, Section 2) for a summary with various pointers to the literature.

However, for general $b$ and $\sigma$, the stationary distribution
$\rho_{\infty}(\cdot)$ does not have a closed-form expression. One good
exception is $b(t,x)=-\nabla f(x)$ and $\sigma(t,x)=\sqrt{2\tau}\,I_{d}$,
where $X$ is governed by the overdamped Langevin equation:

$dX_{t}=-\nabla f(X_{t})dt+\sqrt{2\tau}\,dB_{t},\quad
X_{0}\stackrel{{\scriptstyle d}}{{=}}\mu_{0}(dx).$ (14)

Such a process is time-reversible, and the stationary distribution, under some
growth condition on $f$, is the Gibbs measure

$\nu_{\tau}(dx)=\frac{1}{Z_{\tau}}\exp\left(-\frac{f(x)}{\tau}\right)dx,$ (15)

where $Z_{\tau}:=\int_{\mathbb{R}^{d}}\exp(-f(x)/\tau)dx$ is the normalizing
constant. Much is known about the overdamped Langevin dynamics. For instance,
if $f$ is $\lambda$-convex (i.e. $\nabla^{2}f+\lambda I_{d}$ is positive
definite), the overdamped Langevin process (14) converges exponentially in the
Wasserstein metric with rate $\lambda$ to the Gibbs measure $\nu_{\tau}$
Carrillo et al. (2006). See also Bakry et al. (2014) for modern techniques to
analyze the evolution of the overdamped Langevin equation and generalizations.

Now we turn to the SA process (1). The difference between the overdamped
Langevin process (14) and the process (1) is that the temperature $\tau_{t}$
of the latter is decreasing in time. Due to the time dependence, the limiting
distribution of SA is unknown. As we will see in Section 4, the idea is to
approximate (1) by a process of Gibbs measures with temperature $\tau_{t}$.
Since $\tau_{t}$ decreases to $0$ in the limit, the problem boils down to
studying Gibbs measures at low temperatures. In the next section, we recall
some results of Gibbs measures at low temperatures, which are motivated by
applications in molecular dynamics and Bayesian statistics.

###  3.2. Functional inequalities and the Erying-Kramers law

Here we present functional inequalities of Gibbs measures at low temperatures
$(\tau\to 0)$. Let $\mu$ and $\nu$ be two probability measures on
$\mathbb{R}^{d}$ such that $\mu$ is absolutely continuous relative to $\nu$,
with $d\mu/d\nu$ the Radon-Nikodym derivative. Define the relative entropy or
KL-divergence $H(\mu|\nu)$ of $\mu$ with respect to $\nu$ by

$H(\mu|\nu):=\int\log\bigg{(}\frac{d\mu}{d\nu}\bigg{)}d\mu=\int\frac{d\mu}{d\nu}\log\bigg{(}\frac{d\mu}{d\nu}\bigg{)}d\nu,$
(16)

and the Fisher information $I(\mu|\nu)$ of $\mu$ with respect to $\nu$ by

$I(\mu|\nu):=\frac{1}{2}\int\bigg{|}\nabla\bigg{(}\frac{d\mu}{d\nu}\bigg{)}\bigg{|}^{2}\bigg{(}\frac{d\mu}{d\nu}\bigg{)}^{-1}d\nu.$
(17)

We say that the probability measure $\nu$ satisfies the log-Sobolev inequality
(LSI) with constant $\alpha>0$, if for all probability measures $\mu$ with
$I(\mu|\nu)<\infty$,

$H(\mu|\nu)\leq\frac{1}{\alpha}I(\mu|\nu).$ (18)

The constant $\alpha$ is called the LSI constant for the probability measure
$\nu$. For instance, the LSI constant $\alpha=1$ for $\nu$ the multivariate
Gaussian with mean $0$ and covariance matrix $I_{d}$.

The LSI also plays an important role in studying the convergence rate of the
overdamped Langevin equation. Recall that $\nu_{\tau}$ is the Gibbs measure
defined by (15), and assume that $\nu_{\tau}$ satisfies the LSI with constant
$\alpha_{\tau}>0$. It follows from (Schlichting, 2012, Theorem 1.7) that by
letting $\mu_{\tau,t}$ be the probability distribution of $X_{t}$ defined by
(14), we have

$H(\mu_{\tau,t}|\nu_{\tau})\leq
e^{-2\tau\alpha_{\tau}t}H(\mu_{\tau,0}|\nu_{\tau}).$ (19)

So larger the value of $\alpha_{\tau}$ is, faster the convergence of the
overdamped Langevin process in the KL divergence is. The subscript ‘$\tau$’ in
$\alpha_{\tau}$ suggests the dependence of the LSI constant on the temperature
$\tau$, and we are interested in the asymptotics of $\alpha_{\tau}$ at low
temperatures as $\tau\to 0$. This problem was considered by (Menz and
Schlichting, 2014, Corollary 3.18), who derived a sharp lower bound for
$\alpha_{\tau}$ as $\tau\to 0$.

######  Lemma 1.

Let $f$ satisfy Assumption 1 $\&$ 2. Then the Gibbs measure $\nu_{\tau}$
defined by (15) satisfies the LSI with constant $\alpha_{\tau}>0$ such that

$\alpha_{\tau}\sim C\exp\left(-\frac{E_{*}}{\tau}\right)\quad\mbox{as }\tau\to
0,$ (20)

where $C>0$ depends on $f,d$.

The Eyring-Kramers law provides an estimate on the spectral gap of the
overdamped Langevin equation (14). It dates back to Eyring (1935); Kramers
(1940) in the study of metastability in chemical reactions, and is proved
rigorously by Bovier et al. (2004, 2005). Lemma 1 is the LSI version of the
Eyring-Kramers law, which is stronger than the spectral gap estimate implied
by the Poincaré inequality.

Further define the Wasserstein distance $W_{2}(\mu,\nu)$ between $\mu$ and
$\nu$ by

$W_{2}(\mu,\nu):=\inf_{\Pi}\sqrt{\int|x-y|^{2}\Pi(dx,dy)},$ (21)

where the infimum is over all joint distributions $\Pi$ coupling $\mu$ and
$\nu$. We say that the probability measure $\nu$ satisfies Talagrand’s
inequality with constant $\gamma>0$, if for all probability measure $\mu$ with
$H(\mu|\nu)<\infty$,

$W_{2}(\mu,\nu)\leq\frac{2}{\gamma}H(\mu|\nu).$ (22)

It follows from (Otto and Villani, 2000, Theorem 1) that the LSI implies
Talagrand’s inequality with the same constant. That is, if $\nu$ satisfies the
LSI with constant $\alpha>0$, then $\nu$ also satisfies Talagrand’s inequality
with constant $\gamma=\alpha$. Combining with Lemma 1, we get a lower bound
estimate of Talagrand’s inequality constant for the Gibbs measure
$\nu_{\tau}$.

######  Lemma 2.

Let $f$ satisfy Assumption 1 $\&$ 2. Then the Gibbs measure $\nu_{\tau}$
defined by (15) satisfies Talagrand’s inequality with constant
$\gamma_{\tau}>0$ such that

$\gamma_{\tau}\sim C\exp\left(-\frac{E_{*}}{\tau}\right)\quad\mbox{as }\tau\to
0$ (23)

where $C>0$ depends on $f,d$.

##  4\. Continuous-time simulated annealing

In this section, we prove Theorem 1 by using the ideas developed in Miclo
(1992); Monmarché (2018); Menz et al. (2018). Let $\mu_{t}$ be the probability
measure of $X_{t}$ defined by (1). The key idea is to compare $\mu_{t}$ with
the time-dependent Gibbs measure $\nu_{\tau_{t}}$ given by

$\nu_{\tau_{t}}(dx)=\frac{1}{Z_{\tau_{t}}}\exp\left(-\frac{f(x)}{\tau_{t}}\right)dx,$
(24)

where $Z_{\tau_{t}}:=\int_{\mathbb{R}^{d}}\exp(-f(x)/\tau_{t})$ is the
normalizing constant. Note that $\nu_{\tau_{t}}$ will concentrate on the
minimum point of $f$ as $t\to\infty$ since $\tau_{t}\to 0$ as $t\to\infty$. We
will see that $\nu_{\tau_{t}}$ is close to $\mu_{t}$ in some sense as
$t\to\infty$. The proof of Theorem 1 is broken into four steps.

Step 1: Reduce $\mu_{t}$ to $\nu_{\tau_{t}}$. We establish a bound that
relates $\nu_{\tau_{t}}$ to $\mu_{t}$. Let $(\widetilde{X}_{t};\,t\geq 0)$ be
a process whose distribution is $\nu_{\tau_{t}}$ at time $t$, coupled with
$(X_{t};\,t\geq 0)$ on the same probability space. Fix $\delta>0$. We have

$\displaystyle\mathbb{P}(f(X_{t})>\delta)$
$\displaystyle=\mathbb{P}(f(X_{t})>\delta,f(\widetilde{X}_{t})>\delta)+\mathbb{P}(f(X_{t})>\delta,f(\widetilde{X}_{t})\leq\delta)$
$\displaystyle\leq\mathbb{P}(f(\widetilde{X}_{t})>\delta)+||X_{t}-\widetilde{X}_{t}||_{TV}$
$\displaystyle\leq\mathbb{P}(f(\widetilde{X}_{t})>\delta)+\sqrt{2H(\mu_{t}|\nu_{\tau_{t}})},$
(25)

where we use Pinsker’s inequality (Tsybakov, 2009, Lemma 2.5) in the last
inequation. Now the problem boils down to estimating
$\mathbb{P}(f(\widetilde{X}_{t})>\delta)$ and $H(\mu_{t}|\nu_{\tau_{t}})$.

Step 2: Long-time behavior of $f(\widetilde{X}_{t})$. We study the asymptotics
of $\mathbb{P}(f(\widetilde{X}_{t})>\delta)$ as $t\to\infty$. The following
lemma provides a quantitative estimate of how $\nu_{\tau_{t}}$, or
equivalently $\widetilde{X}_{t}$ concentrates on the minimum point of $f$ as
$t\to\infty$.

######  Lemma 3.

Let $f$ satisfy Assumption 1 $\&$ 2. Assume that $\tau_{t}\sim\frac{E}{\ln t}$
as $t\to\infty$ with $E>E_{*}$. For each $\varepsilon\in(0,\delta)$, there
exist $C>0$ (depending on $\delta,\varepsilon,f,E,d$) such that

$\mathbb{P}(f(\widetilde{X}_{t})>\delta)\leq
Ct^{-\frac{\delta-\varepsilon}{E}}.$ (26)

###### Proof.

Note that

$\mathbb{P}(f(\widetilde{X}_{t})>\delta)=\frac{\int_{f(x)>\delta}\exp(-f(x)/\tau_{t})dx}{\int_{\mathbb{R}^{d}}\exp(-f(x)/\tau_{t})dx}.$
(27)

Under Assumption 1, $f$ has quadratic growth, so at least linear growth at
infinity (Menz and Schlichting, 2014, Lemma 3.14): there exists $C>0$ such
that for $R$ large enough,

$f(x)\geq\min_{|z|=R}f(z)+C(|x|-R)\quad\mbox{for }|x|>R.$

We can also choose $R$ sufficiently large so that $\min_{|z|=R}f(z)>\delta$.
Consequently,

$\displaystyle\int_{f(x)>\delta}\exp(-f(x)/\tau_{t})dx$
$\displaystyle=\int_{f(x)>\delta,|x|\leq
R}\exp(-f(x)/\tau_{t})dx+\int_{f(x)>\delta,|x|>R}\exp(-f(x)/\tau_{t})dx$
$\displaystyle=e^{-\frac{\delta}{\tau_{t}}}\operatorname{Vol}(B_{R})(1+\mathcal{O}(\tau_{t})),$
(28)

where $\operatorname{Vol}(B_{R})$ is the volume of a ball with radius $R$.
Further by Laplace’s method,

$\int_{\mathbb{R}^{d}}\exp(-f(x)/\tau_{t})dx\sim C(\tau_{t})^{\frac{d}{2}}.$
(29)

By injecting (4), (29) into (27),we get

$\mathbb{P}(f(\widetilde{X}_{t})>\delta)\leq Ct^{-\frac{\delta}{E}}(\ln
t)^{\frac{d}{2}},$ (30)

which clearly yields (26) since $(\ln
t)^{\frac{d}{2}}/t^{\frac{\varepsilon}{E}}\to 0$ as $t\to\infty$. ∎

######  Remark 1.

It is interesting to get a bound for $\mathbb{P}(f(\widetilde{X}_{t})>\delta)$
when the dimension $d$ is large. As mentioned in the introduction, the Laplace
bound (29) may fail when $d,t\to\infty$ simultaneously. Recall that $m_{0}$ is
the minimum point of $f$. By continuity of $f$, there exists $r>0$ such that
$f(x)<\varepsilon$ when $|x-m_{0}|<r$. Thus,

$\displaystyle\int_{\mathbb{R}^{d}}\exp(-f(x)/\tau_{t})dx$
$\displaystyle\geq\int_{|x-m_{0}|<r}\exp(-f(x)/\tau_{t})dx$ $\displaystyle\geq
e^{-\frac{\varepsilon}{\tau_{t}}}\operatorname{Vol}(B_{r}).$ (31)

Further, if $t/e^{\frac{Ed}{CR}}\to\infty$ as $d\to\infty$,

$\int_{f(x)>\delta}\exp(-f(x)/\tau_{t})dx=e^{-\frac{\delta}{\tau_{t}}}\operatorname{Vol}(B_{R})(1+\mathcal{O}(\tau_{t}d)),$
(32)

Combining (1) and (32), we get

$\mathbb{P}(f(\widetilde{X}_{t})>\delta)\leq
C\gamma^{d}t^{-\frac{\delta-\varepsilon}{E}},$ (33)

where $C>0$ depends on $\delta,\varepsilon,f,E$, and
$\gamma=\max(R/r,e^{\frac{\delta-\varepsilon}{CR}})$. Also note that Raginsky
et al. (2017) obtained the bound $\mathbb{E}f(\widetilde{X}_{t})\leq Cd/\ln
t$. By Markov’s inequality, we get

$\mathbb{P}(f(\widetilde{X}_{t})>\delta)\leq C\delta^{-1}d\,(\ln t)^{-1}.$
(34)

In comparison with (34), the bound (33) is better in ‘$t$’ but worse in ‘$d$’.
In terms of relaxation time, i.e. letting
$\mathbb{P}(f(\widetilde{X}_{t})>\delta)$ be of constant order, both estimates
show an exponential dependence of $t$ on $d$. This suggests that SA is
exponentially slow as the dimension increases.

Step 3: Differential inequality for $H(\mu_{t}|\nu_{\tau_{t}})$. To get an
estimate of $H(\mu_{t}|\nu_{\tau_{t}})$, we need to consider the time
derivative $\frac{d}{dt}H(\mu_{t}|\nu_{\tau_{t}})$. The following lemma is a
reformulation of (Miclo, 1992, Proposition 3). For ease of reference, we give
a simplified proof here. First let us convent some notation. For an absolutely
continuous measure $\mu(dx)$, we abuse the notation $\mu(dx)=\mu(x)dx$, i.e.
$\mu(x)$ is the density of $\mu(dx)$. So for two such probability measures
$\mu$ and $\nu$, the Radon-Nikodym derivative $\frac{d\mu}{d\nu}(x)$ is
identified with $\frac{\mu(x)}{\nu(x)}$.

######  Lemma 4.

Let $\tau_{t}$ be decreasing in $t$. We have

$\frac{d}{dt}H(\mu_{t}|\nu_{\tau_{t}})\leq-2\tau_{t}I\left(\mu_{t}|\nu_{\tau_{t}}\right)+\frac{d}{dt}\bigg{(}\frac{1}{\tau_{t}}\bigg{)}\,\mathbb{E}f(X_{t}).$
(35)

where $I(\mu_{t}|\nu_{\tau_{t}})$ is the Fisher information defined by (17).

###### Proof.

Observe that

$\displaystyle\frac{d}{dt}H(\mu_{t}|\nu_{\tau_{t}})$
$\displaystyle=\frac{d}{dt}\int\mu_{t}\ln\bigg{(}\frac{\mu_{t}}{\nu_{\tau_{t}}}\bigg{)}dx$
$\displaystyle=\underbrace{\int\frac{\partial\mu_{t}}{\partial
t}\ln\bigg{(}\frac{\mu_{t}}{\nu_{\tau_{t}}}\bigg{)}dx}_{(a)}+\underbrace{\int\mu_{t}\frac{\partial}{\partial
t}\bigg{(}\ln\bigg{(}\frac{\mu_{t}}{\nu_{\tau_{t}}}\bigg{)}\bigg{)}dx}_{(b)}.$
(36)

We first consider the term (a). Recall that $\mu_{t}$ satisfies the Fokker-
Planck equation (13). Together with the fact that
$\nabla(\tau_{t}\nu_{\tau_{t}})=-\nu_{\tau_{t}}\nabla f$, we have

$\frac{\partial\mu_{t}}{\partial
t}=\nabla\cdot\bigg{(}\tau_{t}\nu_{\tau_{t}}\nabla\bigg{(}\frac{\mu_{t}}{\nu_{\tau_{t}}}\bigg{)}\bigg{)}.$
(37)

By injecting (37) into the term (a) and further by integration by parts, we
get

$\displaystyle(a)$
$\displaystyle=\int\nabla\cdot\bigg{(}\tau_{t}\nu_{\tau_{t}}\nabla\bigg{(}\frac{\mu_{t}}{\nu_{\tau_{t}}}\bigg{)}\bigg{)}\ln\bigg{(}\frac{\mu_{t}}{\nu_{\tau_{t}}}\bigg{)}dx$
$\displaystyle=-\int\tau_{t}\nu_{\tau_{t}}\nabla\bigg{(}\frac{\mu_{t}}{\nu_{\tau_{t}}}\bigg{)}\cdot\nabla\ln\bigg{(}\frac{\mu_{t}}{\nu_{\tau_{t}}}\bigg{)}dx$
$\displaystyle=-\tau_{t}\int\bigg{|}\nabla\bigg{(}\frac{\mu_{t}}{\nu_{\tau_{t}}}\bigg{)}\bigg{|}^{2}\bigg{(}\frac{\mu_{t}}{\nu_{\tau_{t}}}\bigg{)}^{-1}d\nu_{\tau_{t}}=-2\tau_{t}I(\mu_{t}|\nu_{\tau_{t}}).$
(38)

Now we consider the term (b). Direct computation leads to

$\displaystyle(b)=\int\bigg{(}\frac{\partial\mu_{t}}{\partial
t}-\frac{\mu_{t}}{\nu_{\tau_{t}}}\frac{\partial\nu_{\tau_{t}}}{\partial
t}\bigg{)}dx$ $\displaystyle=-\int\frac{\partial}{\partial
t}\left(\ln\nu_{\tau_{t}}\right)d\mu_{t}$
$\displaystyle=\int\frac{d}{dt}\left(\ln
Z_{\tau_{t}}\right)d\mu_{t}+\frac{d}{dt}\bigg{(}\frac{1}{\tau_{t}}\bigg{)}\mathbb{E}f(X_{t})$
$\displaystyle\leq\frac{d}{dt}\bigg{(}\frac{1}{\tau_{t}}\bigg{)}\,\mathbb{E}f(X_{t}),$
(39)

where we use the fact that $\int\mu_{t}dx=1$ in the second equation, and that
$\tau_{t}$ is decreasing in $t$ in the last inequality. Combining (4) with (4)
and (4) yields (35). ∎

Step 4: Estimating $H(\mu_{t}|\nu_{\tau_{t}})$ via the Erying-Kramers law.
Note that there are two terms on the right hand side of (35). We start with an
estimate of the second term.

######  Lemma 5.

Let $f$ satisfy Assumption 1, and assume that the condition (4) for $\mu_{0}$
holds. For each $\varepsilon>0$, there exists $C>0$ (depending on
$\varepsilon,f$) such that

$\mathbb{E}f(X_{t})\leq C(1+t)^{\varepsilon}.$ (40)

###### Proof.

It is easy to see that Assumption 1 implies Assumption $\mbox{H}_{1}$ in Miclo
(1992). Together with the moment condition (4), the proof follows the line of
reasoning in (Miclo, 1992, Lemma 2). ∎

Now we apply the Eyring-Kramers law, combining with a Grönwall-type argument
to bound $H(\mu_{t}|\nu_{\tau_{t}})$ for large $t$.

######  Lemma 6.

Let $f$ satisfy Assumption 1 $\&$ 2, and assume that the condition (4) for
$\mu_{0}$ holds. Assume that $\tau_{t}\sim\frac{E}{\ln t}$ and
$\frac{d}{dt}\left(\frac{1}{\tau_{t}}\right)=\mathcal{O}\left(\frac{1}{t}\right)$
as $t\to\infty$ with $E>E_{*}$. For each $\varepsilon>0$, there exists $C>0$
(depending on $\varepsilon,f,E,d$) such that

$H(\mu_{t}|\nu_{\tau_{t}})\leq
Ct^{-\left(1-\frac{E_{*}}{E}-\varepsilon\right)}.$ (41)

###### Proof.

Using Lemma 1 and the bound (35), we have

$\frac{d}{dt}H(\mu_{t}|\nu_{\tau_{t}})\leq-2\tau_{t}\alpha_{t}H(\mu_{t}|\nu_{\tau_{t}})+\frac{C}{t}\,\mathbb{E}f(X_{t}),$
(42)

where $\alpha_{t}$ is the LSI constant for the Gibbs measure $\nu_{\tau_{t}}$.
By the Eyring-Kramers formula (20), for each $\varepsilon>0$, there exist
$C>0$ and $t_{0}>0$,

$2\tau_{t}\alpha_{t}\geq
Ct^{-\left(\frac{E_{*}}{E}-\varepsilon\right)}\quad\mbox{for }t\geq t_{0}.$
(43)

Combining (42) with (40), (43), we get

$\frac{d}{dt}H(\mu_{t}|\nu_{\tau_{t}})\leq-
Ct^{-\left(\frac{E_{*}}{E}-\varepsilon\right)}H(\mu_{t}|\nu_{\tau_{t}})+C^{\prime}t^{-1+\varepsilon}.$
(44)

Fix $\varepsilon\in(0,\frac{1}{2}-\frac{E_{*}}{2E})$, let

$Q(t):=H(\mu_{t}|\nu_{\tau_{t}})-\frac{2C^{\prime}}{C}t^{-1+\frac{E_{*}}{E}+2\varepsilon}.$

Then for $t_{0}$ sufficiently large and $t\geq t_{0}$, we have
$\frac{d}{dt}Q(t)\leq-Ct^{-\frac{E_{*}}{E}+\varepsilon}Q(t)$ by (44). This
implies that $Q(t)\leq
Q(t_{0})e^{-C\int_{t_{0}}^{t}s^{-\frac{E_{*}}{E}+\varepsilon}ds}$. Thus,

$H(\mu_{t}|\nu_{\tau_{t}})\leq\frac{2C^{\prime}}{C}t^{-1+\frac{E_{*}}{E}+2\varepsilon}+H(\mu_{t_{0}}|\nu_{t_{0}})e^{-\frac{C}{\kappa}(t^{\kappa}-t_{0}^{\kappa})},$
(45)

where $\kappa:=1-\frac{E_{*}}{E}-\varepsilon>0$. Note that the first term on
the right hand side of (45) dominates, and the conclusion follows. ∎

Finally, by injecting (26), (41) into (4), we get the desired estimate (5).

##  5\. Discrete simulated annealing

This section is devoted to the proof of Theorem 2. The idea is close to that
employed for the continuous-time SA process (1). However, the analysis is more
complicated due to discretization, and additional tools from Vempala and
Wibisono (2019) on the convergence of ULA are used. Recall that $\eta_{k}$ is
the step size at iteration $k$, and $\Theta_{k}:=\sum_{j\leq k}\eta_{j}$ is
the cumulative step size up to iteration $k$. Let $\mu_{k}$ be the probability
density of $x_{k}$ defined by (2), and

$\nu_{\tau_{\scaleto{\Theta_{k}}{4pt}}}(dx)=\frac{1}{Z_{\tau_{\scaleto{\Theta_{k}}{4pt}}}}\exp\left(-\frac{f(x)}{\tau_{\scaleto{\Theta_{k}}{5pt}}}\right)dx,$
(46)

where
$Z_{\tau_{\scaleto{\Theta_{k}}{4pt}}}:=\int_{\mathbb{R}^{d}}\exp(-f(x)/\tau_{\scaleto{\Theta_{k}}{5pt}})dx$
is the normalizing constant. We divide the proof into four steps.

Step 1: Reduce $\mu_{k}$ to $\nu_{\tau_{\scaleto{\Theta_{k}}{4pt}}}$. This
step is similar to Step $1$ $\&$ $2$ in Section 4. Let
$(\widetilde{x}_{k};\,k\geq 0)$ be a sequence whose distribution is
$\nu_{\tau_{\scaleto{\Theta_{k}}{4pt}}}$ at epoch $k$, coupled with
$(x_{k};\,k\geq 0)$ on the same probability space. Fix $\delta>0$. The same
argument as in (4) shows that

$\mathbb{P}(f(x_{k})>\delta)\leq\mathbb{P}(f(\widetilde{x}_{k})>\delta)+\sqrt{2H(\mu_{k}|\nu_{\tau_{\scaleto{\Theta_{k}}{4pt}}})}.$
(47)

Assume that $\Theta_{k}\to\infty$ and
$\tau_{\scaleto{\Theta_{k}}{5pt}}\sim\frac{E}{\ln\Theta_{k}}$ as $k\to\infty$
with $E>E_{*}$. By Lemma 3, we get a bound for the first term on the right
hand side of (47). That is, for each $\varepsilon\in(0,\delta)$, there exists
$C>0$ (depending on $\varepsilon,\delta,f,E,d$) such that

$\mathbb{P}(f(\widetilde{x}_{k})>\delta)\leq
C\Theta_{k}^{-\frac{\delta-\varepsilon}{E}}.$ (48)

So it remains to estimate $H(\mu_{k}|\nu_{\tau_{\scaleto{\Theta_{k}}{4pt}}})$,
which is the task of the next three steps.

Step 2: Continuous-time coupling. To make use of continuous-time tools, we
couple the sequence $(x_{k};\,k\geq 0)$ by a continuous-time process
$(X_{t};\,t\geq 0)$ such that $(X_{\Theta_{k}};\,k\geq 0)$ has the same
distribution as $(x_{k};\,k\geq 0)$. To do this, define the process $X$ by

$dX_{t}=-\nabla
f(x_{k})dt+\sqrt{2\tau_{\scaleto{\Theta_{k}}{5pt}}}dB_{t},\quad
t\in[\Theta_{k},\Theta_{k+1}),$ (49)

where we identify $X_{\Theta_{k}}$ with $x_{k}$. Note that the Fokker-Planck
equation (13) plays an important role in the analysis of continuous-time SA
(1). It is desirable to get a version of the Fokker-Planck equation for the
coupled process (49). The result is stated as follows.

######  Lemma 7.

For $t\in[\Theta_{k},\Theta_{k+1})$, the probability density $\mu_{t}$ of
$X_{t}$ defined by (49) satisfies the following equation:

$\frac{\partial\mu_{t}}{\partial
t}=\nabla\cdot\bigg{(}\tau_{\scaleto{\Theta_{k}}{5pt}}\nu_{\tau_{\scaleto{\Theta_{k}}{4pt}}}\nabla\bigg{(}\frac{\mu_{t}}{\nu_{\tau_{\scaleto{\Theta_{k}}{4pt}}}}\bigg{)}\bigg{)}+\nabla\cdot\left(\mu_{t}\,\mathbb{E}[\nabla
f(x_{k})-\nabla f(X_{t})|X_{t}=x]\right).$ (50)

###### Proof.

Let $\mu_{t|s}(x|y)$ be the conditional probability
$\mathbb{P}(X_{t}=x|X_{s}=y)$. By conditioning on $X_{\Theta_{k}}=x_{k}$, we
have

$\frac{\partial\mu_{t|\Theta_{k}}(x|x_{k})}{\partial
t}=\nabla\cdot(\mu_{t|\Theta_{k}}(x|x_{k})\nabla
f(x_{k}))+\tau_{\Theta_{k}}\Delta\mu_{t|\Theta_{k}}(x|x_{k}).$ (51)

By integrating (51) against $\mu_{\Theta_{k}}$, and using the fact that
$\mu_{t|\Theta_{k}}(x|x_{k})\mu_{\Theta_{k}}(x_{k})=\mu_{t}(x)\mu_{\Theta_{k}|t}(x_{k}|x)$,
we get

$\frac{\partial\mu_{t}}{\partial t}=\nabla\cdot(\mu_{t}(x)\,\mathbb{E}[\nabla
f(x_{k})|X_{t}=x])+\tau_{\scaleto{\Theta_{k}}{5pt}}\Delta\mu_{t}.$ (52)

Further by (37), we have

$\nabla\cdot\bigg{(}\tau_{\scaleto{\Theta_{k}}{5pt}}\nu_{\tau_{\scaleto{\Theta_{k}}{4pt}}}\nabla\bigg{(}\frac{\mu_{t}}{\nu_{\tau_{\scaleto{\Theta_{k}}{4pt}}}}\bigg{)}\bigg{)}=\nabla\cdot(\mu_{t}\nabla
f(x))+\tau_{\scaleto{\Theta_{k}}{5pt}}\Delta\mu_{t}.$ (53)

Combining (52) and (53) yields (50). ∎

There are two terms on the right hand side of (50). Comparing to (37), the
first term is the usual Fokker-Planck term, while the second term corresponds
to the discretization error.

Step 3: One-step analysis of
$H(\mu_{k}|\nu_{\tau_{\scaleto{\Theta_{k}}{4pt}}})$. Here we use the coupled
process (49) to study the one-step decay of
$H(\mu_{k}|\nu_{\tau_{\scaleto{\Theta_{k}}{4pt}}})$.

######  Lemma 8.

Let $f$ satisfy Assumption 1, 2 $\&$ 3, and assume that the condition (4) for
$\mu_{0}$ holds. Assume that $\tau_{t}\sim\frac{E}{\ln t}$ and
$\frac{d}{dt}\left(\frac{1}{\tau_{t}}\right)=\mathcal{O}\left(\frac{1}{t}\right)$
as $t\to\infty$ with $E>E_{*}$. Also assume that $\Theta_{k}\to\infty$ and
$\eta_{k+1}\Theta_{k}\to 0$ as $k\to\infty$. Then, for each $\varepsilon>0$,
there exist $C,C^{\prime}>0$ (depending on $\varepsilon,f,E,d$) such that

$\displaystyle H(\mu_{k+1}|\nu_{\tau_{\scaleto{\Theta_{k+1}}{4pt}}})$
$\displaystyle\leq\bigg{(}1-C\eta_{k+1}\Theta_{k}^{-(\frac{E_{*}}{E}-\varepsilon)}\bigg{)}H(\mu_{k}|\nu_{\tau_{\scaleto{\Theta_{k}}{4pt}}})$
$\displaystyle\qquad\qquad\qquad\quad+C^{\prime}(\eta_{k+1}^{2}+\eta_{k+1}^{3}\ln\Theta_{k}+\eta_{k+1}\Theta_{k}^{-1+\varepsilon}).$
(54)

###### Proof.

Write

$H(\mu_{k+1}|\nu_{\tau_{\scaleto{\Theta_{k+1}}{4pt}}})=\underbrace{H(\mu_{k+1}|\nu_{\tau_{\scaleto{\Theta_{k}}{4pt}}})}_{(a)}+\underbrace{(H(\mu_{k+1}|\nu_{\tau_{\scaleto{\Theta_{k+1}}{4pt}}})-H(\mu_{k+1}|\nu_{\tau_{\scaleto{\Theta_{k}}{4pt}}}))}_{(b)}.$
(55)

We first use the coupled process (49) to study the term $(a)$. Note that

$\displaystyle\frac{d}{dt}H(\mu_{t}|\nu_{\tau_{\scaleto{\Theta_{k}}{4pt}}})$
$\displaystyle=\int\frac{\partial\mu_{t}}{\partial
t}\ln\bigg{(}\frac{\mu_{t}}{\nu_{\tau_{\scaleto{\Theta_{k}}{4pt}}}}\bigg{)}dx+\int\mu_{t}\frac{d}{dt}\ln\bigg{(}\frac{\mu_{t}}{\nu_{\tau_{\scaleto{\Theta_{k}}{4pt}}}}\bigg{)}dx$
$\displaystyle=\int\nabla\cdot\bigg{(}\tau_{\scaleto{\Theta_{k}}{5pt}}\nu_{\tau_{\scaleto{\Theta_{k}}{4pt}}}\nabla\bigg{(}\frac{\mu_{t}}{\nu_{\tau_{\scaleto{\Theta_{k}}{4pt}}}}\bigg{)}\bigg{)}\ln\bigg{(}\frac{\mu_{t}}{\nu_{\tau_{\scaleto{\Theta_{k}}{4pt}}}}\bigg{)}dx$
$\displaystyle\quad+\underbrace{\int\nabla\cdot\left(\mu_{t}\,\mathbb{E}[\nabla
f(x_{k})-\nabla
f(X_{t})|X_{t}=x]\right)\ln\bigg{(}\frac{\mu_{t}}{\nu_{\tau_{\scaleto{\Theta_{k}}{4pt}}}}\bigg{)}dx}_{(c)}+\frac{d}{dt}\int\mu_{t}(dx)$
$\displaystyle=-2\tau_{\scaleto{\Theta_{k}}{5pt}}I(\mu_{t}|\nu_{\tau_{\scaleto{\Theta_{k}}{4pt}}})+(c),$
(56)

where we use (50) in the second equation, and (4) in the third equation. Now
we need to estimate the term $(c)$ in (5). By integration by parts and the
fact that $a\cdot
b\leq\frac{1}{\tau_{\scaleto{\Theta_{k}}{4pt}}}|a|^{2}+\frac{\tau_{\scaleto{\Theta_{k}}{4pt}}}{4}|b|^{2}$,
we get

$\displaystyle(c)$ $\displaystyle=\mathbb{E}\bigg{(}(\nabla f(X_{t})-\nabla
f(x_{k}))\cdot\nabla\log\bigg{(}\frac{\mu_{t}}{\nu_{\tau_{\scaleto{\Theta_{k}}{4pt}}}}\bigg{)}\bigg{)}$
$\displaystyle\leq\frac{1}{\tau_{\scaleto{\Theta_{k}}{5pt}}}\mathbb{E}|\nabla
f(X_{t})-\nabla
f(x_{k}))|^{2}+\frac{\tau_{\scaleto{\Theta_{k}}{5pt}}}{4}\mathbb{E}\bigg{|}\log\bigg{(}\frac{\mu_{t}}{\nu_{\tau_{\scaleto{\Theta_{k}}{4pt}}}}\bigg{)}\bigg{|}^{2}$
$\displaystyle\leq\frac{L^{2}}{\tau_{\scaleto{\Theta_{k}}{5pt}}}\mathbb{E}|X_{t}-x_{k}|^{2}+\frac{\tau_{\scaleto{\Theta_{k}}{5pt}}}{2}I(\mu_{t}|\nu_{\tau_{\scaleto{\Theta_{k}}{4pt}}}),$
(57)

where $L$ is the Lipschitz constant of $\nabla f$ by Assumption 3. Recall from
(49) that $X_{t}-x_{k}=-\nabla
f(x_{k})(t-\Theta_{k})+\sqrt{2\tau_{\scaleto{\Theta_{k}}{5pt}}(t-\Theta_{k})}Z$,
where $Z$ is standard normal. Consequently,

$\displaystyle\mathbb{E}|X_{t}-x_{k}|^{2}$
$\displaystyle=(t-x_{k})^{2}\mathbb{E}|\nabla
f(x_{k})|^{2}+2\tau_{\scaleto{\Theta_{k}}{5pt}}(t-x_{k})d$
$\displaystyle\leq\eta_{k+1}^{2}\mathbb{E}|\nabla
f(x_{k})|^{2}+C\tau_{\scaleto{\Theta_{k}}{5pt}}\eta_{k+1}.$ (58)

According to Lemma 2, $\nu_{\tau_{\scaleto{\Theta_{k}}{4pt}}}$ satisfies
Talagrand’s inequality with constant
$\gamma_{\tau_{\scaleto{\Theta_{k}}{4pt}}}\sim\kappa\exp(-E_{*}/\tau_{\scaleto{\Theta_{k}}{5pt}})$.
Further by (Vempala and Wibisono, 2019, Lemma 10),

$\mathbb{E}|\nabla
f(x_{k})|^{2}\leq\frac{C}{\gamma_{\tau_{\scaleto{\Theta_{k}}{4pt}}}}H(\mu_{k}|\nu_{\tau_{\scaleto{\Theta_{k}}{4pt}}})+C.$
(59)

Combining (5) with (5), (59) and the fact that
$\tau_{\scaleto{\Theta_{k}}{5pt}}\sim\frac{E}{\ln\Theta_{k}}$ as $k\to\infty$,
we have

$(c)\leq
C\bigg{(}\eta_{k+1}^{2}\Theta_{k}^{\frac{E_{*}}{E}}\ln\Theta_{k}\bigg{)}H(\mu_{k}|\nu_{\tau_{\scaleto{\Theta_{k}}{4pt}}})+C(\eta_{k+1}+\eta_{k+1}^{2}\ln\Theta_{k})+\frac{\tau_{\scaleto{\Theta_{k}}{5pt}}}{2}I(\mu_{t}|\nu_{\tau_{\scaleto{\Theta_{k}}{4pt}}}).$
(60)

Injecting (60) into (5) and further by Lemma 1, we get

$\displaystyle\frac{d}{dt}H(\mu_{t}|\nu_{\tau_{\scaleto{\Theta_{k}}{4pt}}})$
$\displaystyle\leq-\frac{3}{2}\tau_{\scaleto{\Theta_{k}}{5pt}}I(\mu_{t}|\nu_{\tau_{\scaleto{\Theta_{k}}{4pt}}})+C^{\prime}\bigg{(}\eta_{k+1}^{2}\Theta_{k}^{\frac{E_{*}}{E}}\ln\Theta_{k}\bigg{)}H(\mu_{k}|\nu_{\tau_{\scaleto{\Theta_{k}}{4pt}}})+C^{\prime}(\eta_{k+1}+\eta_{k+1}^{2}\ln\Theta_{k})$
$\displaystyle\leq-\frac{3}{2}C\Theta_{k}^{-(\frac{E_{*}}{E}-\varepsilon)}H(\mu_{t}|\nu_{\tau_{\scaleto{\Theta_{k}}{4pt}}})+C^{\prime}\bigg{(}\eta_{k+1}^{2}\Theta_{k}^{\frac{E_{*}}{E}+\varepsilon}H(\mu_{k}|\nu_{\tau_{\scaleto{\Theta_{k}}{4pt}}})+(\eta_{k+1}+\eta_{k+1}^{2}\ln\Theta_{k})\bigg{)},$

Now by a Grönwall argument, we have

$\displaystyle H(\mu_{k+1}|\nu_{\tau_{\scaleto{\Theta_{k}}{4pt}}})$
$\displaystyle\leq
e^{-\frac{3}{2}C\eta_{k+1}\Theta_{k}^{-(\frac{E_{*}}{E}-\varepsilon)}}\left((1+C^{\prime}\eta_{k+1}^{3}\Theta^{\frac{E_{*}}{E}+\varepsilon})H(\mu_{k}|\nu_{\tau_{\scaleto{\Theta_{k}}{4pt}}})+C^{\prime}(\eta_{k+1}^{2}+\eta_{k+1}^{3}\ln\Theta_{k})\right)$
$\displaystyle\leq
e^{-\frac{5}{4}C\eta_{k+1}\Theta_{k}^{-(\frac{E_{*}}{E}-\varepsilon)}}H(\mu_{k}|\nu_{\tau_{\scaleto{\Theta_{k}}{4pt}}})+C^{\prime}(\eta_{k+1}^{2}+\eta_{k+1}^{3}\ln\Theta_{k})$
$\displaystyle\leq\bigg{(}1-C\eta_{k+1}\Theta_{k}^{-(\frac{E_{*}}{E}-\varepsilon)}\bigg{)}H(\mu_{k}|\nu_{\tau_{\scaleto{\Theta_{k}}{4pt}}})+C^{\prime}(\eta_{k+1}^{2}+\eta_{k+1}^{3}\ln\Theta_{k}),$
(61)

where we use the fact that $\eta_{k+1}\Theta_{k}^{\frac{E_{*}}{E}}\to 0$ as
$k\to\infty$ in the second inequality.

Now we consider the term $(b)$ in (55). Note that

$\displaystyle
H(\mu_{k+1}|\nu_{\tau_{\scaleto{\Theta_{k+1}}{4pt}}})-H(\mu_{k+1}|\nu_{\tau_{\scaleto{\Theta_{k}}{4pt}}})$
$\displaystyle=\ln\bigg{(}\frac{Z_{\tau_{\scaleto{\Theta_{k+1}}{4pt}}}}{Z_{\tau_{\scaleto{\Theta_{k}}{4pt}}}}\bigg{)}+\bigg{(}\frac{1}{\tau_{\scaleto{\Theta_{k+1}}{5pt}}}-\frac{1}{\tau_{\scaleto{\Theta_{k}}{4pt}}}\bigg{)}\mathbb{E}f(x_{k+1})$
$\displaystyle\leq C\frac{\eta_{k+1}}{\Theta_{k}}\mathbb{E}f(x_{k+1}).$ (62)

We claim that for each $\varepsilon>0$, $\mathbb{E}f(x_{k+1})\leq C\leq
C\Theta_{k}^{\varepsilon}$. Choose $C>0$ sufficiently large, and let
$\mathbb{E}f(x_{k+1})$ be the first term exceeding $C$. By Assumption 3,

$\displaystyle f(x_{k+1})\leq f(x_{k})-\eta_{k}|\nabla
f(x_{k})|^{2}+\sqrt{2\tau_{\scaleto{\Theta_{k}}{5pt}}\eta_{k}}\nabla
f(x_{k})\cdot Z_{k}+\frac{L}{2}|\eta_{k}\nabla
f(x_{k})+\sqrt{2\tau_{\scaleto{\Theta_{k}}{5pt}}\eta_{k}}Z_{k}|^{2}$

Further by taking expectation, we get

$\mathbb{E}f(x_{k+1})-\mathbb{E}f(x_{k})\leq-\eta_{k}\bigg{(}1-\frac{\eta_{k}L}{2}\bigg{)}\mathbb{E}|\nabla
f(x_{k})|^{2}+Ld\tau_{\scaleto{\Theta_{k}}{5pt}}\eta_{k}.$ (63)

Thus, $\mathbb{E}f(x_{k+1})-\mathbb{E}f(x_{k})\leq
Ld\tau_{\scaleto{\Theta_{k}}{5pt}}\eta_{k}$ which implies that
$\mathbb{E}f(x_{k})>C-1$ for $k$ large enough. By Assumption 1(ii),
$\mathbb{E}|\nabla f(x_{k})|^{2}>C^{\prime}$ for some $C^{\prime}>0$.
Combining with (63), we have $\mathbb{E}f(x_{k})>\mathbb{E}f(x_{k+1})\geq C$
for $k$ large enough. This contradicts the fact that $\mathbb{E}f(x_{k+1})$ is
the first term exceeding $C$. Now by (5), we get

$H(\mu_{k+1}|\nu_{\tau_{\scaleto{\Theta_{k+1}}{4pt}}})-H(\mu_{k+1}|\nu_{\tau_{\scaleto{\Theta_{k}}{4pt}}})\leq
C\eta_{k+1}\Theta_{k}^{-1+\varepsilon}.$ (64)

Combining (55) with (5), (64) yields (8). ∎

Step 4: Estimating $H(\mu_{k}|\nu_{\tau_{\scaleto{\Theta_{k}}{4pt}}})$. We use
Lemma 8 to derive an estimate for
$H(\mu_{k}|\nu_{\tau_{\scaleto{\Theta_{k}}{4pt}}})$. Under the condition (6),
the term $\eta_{k+1}\Theta_{k}^{-1+\varepsilon}$ dominates $\eta_{k+1}^{2}$,
$\eta_{k+1}^{3}\ln\Theta_{k}$ as $k\to\infty$. Thus, the recursion (8) yields

$H(\mu_{k+1}|\nu_{\tau_{\scaleto{\Theta_{k+1}}{4pt}}})\leq\bigg{(}1-C\eta_{k+1}\Theta_{k}^{-(\frac{E_{*}}{E}-\varepsilon)}\bigg{)}H(\mu_{k}|\nu_{\tau_{\scaleto{\Theta_{k}}{4pt}}})+C^{\prime}\eta_{k+1}\Theta_{k}^{-1+\varepsilon}.$

Since $E_{*}/E<1$, a similar argument as in Lemma 6 shows that

$H(\mu_{k+1}|\nu_{\tau_{\scaleto{\Theta_{k+1}}{4pt}}})-C\Theta_{k}^{-(1-\frac{E_{*}}{E}-\varepsilon)}\leq\bigg{(}1-C^{\prime}\eta_{k+1}\Theta_{k}^{-(\frac{E_{*}}{E}-\varepsilon)}\bigg{)}\bigg{(}H(\mu_{k}|\nu_{\tau_{\scaleto{\Theta_{k}}{4pt}}})-C\Theta_{k-1}^{-(1-\frac{E_{*}}{E}-\varepsilon)}\bigg{)}.$

This together with the condition (7) implies that

$H(\mu_{k}|\nu_{\tau_{\scaleto{\Theta_{k}}{4pt}}})\leq
C\Theta_{k}^{-(1-\frac{E_{*}}{E}-\varepsilon)}.$ (65)

By injecting (48) and (65) into (47) we obtain (8).

##  6\. Conclusion

In this paper, we study the convergence rate of SA in both continuous and
discrete settings. The main tool is functional inequalities for the Gibbs
measure at low temperatures. We prove that the tail probability, in both
settings, exhibits a polynomial decay in time. The decay rate is also given as
function of the model parameters. In the discrete setting, we derive a
condition on the step size to ensure the convergence to the global minimum.
This condition may be useful in tuning the step size.

There are a few directions to extend this work. For instance, one can study
the convergence rate of SA for Lévy flight with a suitable cooling schedule.
Another problem is to study the dependence of the convergence rate in the
dimension $d$. This requires a deep understanding of the Eyring-Kramers law in
high dimension, and is related to the Laplace approximation of high
dimensional integrals. Both problems are worth exploring, but may be
challenging.

Acknowledgement: Tang thanks Georg Menz for helpful discussions. He also
gratefully acknowledges financial support through a start-up grant at Columbia
University. Zhou gratefully acknowledges financial supports through a start-up
grant at Columbia University and through the Nie Center for Intelligent Asset
Management.

## References

  * Bakry et al. (2014) D. Bakry, I. Gentil, and M. Ledoux.  _Analysis and geometry of Markov diffusion operators_ , volume 348 of _Grundlehren der Mathematischen Wissenschaften_.  Springer, 2014.
  * Bovier et al. (2004) A. Bovier, M. Eckhoff, V. Gayrard, and M. Klein.  Metastability in reversible diffusion processes. I. Sharp asymptotics for capacities and exit times.  _J. Eur. Math. Soc._ , 6(4):399–424, 2004.
  * Bovier et al. (2005) A. Bovier, V. Gayrard, and M. Klein.  Metastability in reversible diffusion processes. II. Precise asymptotics for small eigenvalues.  _J. Eur. Math. Soc._ , 7(1):69–99, 2005.
  * Carrillo et al. (2006) J. A. Carrillo, R. J. McCann, and C. Villani.  Contractions in the 2-Wasserstein length space and thermalization of granular media.  _Arch. Ration. Mech. Anal._ , 179(2):217–263, 2006.
  * Cerny (1985) V. Cerny.  Thermodynamical approach to the traveling salesman problem: an efficient simulation algorithm.  _J. Optim. Theory Appl._ , 45(1):41–51, 1985\.
  * Chen et al. (2020) X. Chen, S. S. Du, and X. T. Tong.  On stationary-point hitting time and ergodicity of stochastic gradient Langevin dynamics.  _J. Mach. Learn. Res._ , 21:Paper No. 68, 41, 2020.
  * Cherny and Engelbert (2005) A. S. Cherny and H.-J. Engelbert.  _Singular stochastic differential equations_ , volume 1858 of _Lecture Notes in Mathematics_.  Springer-Verlag, 2005.
  * Chiang et al. (1987) T.-S. Chiang, C.-R. Hwang, and S. J. Sheu.  Diffusion for global optimization in ${\bf R}^{n}$.  _SIAM J. Control Optim._ , 25(3):737–753, 1987\.
  * Delahaye et al. (2019) D. Delahaye, S. Chaimatanan, and M. Mongeau.  Simulated annealing: from basics to applications.  In _Handbook of metaheuristics_ , volume 272 of _Internat. Ser. Oper. Res. Management Sci._ , pages 1–35. Springer, 2019.
  * Durmus and Moulines (2017) A. Durmus and E. Moulines.  Nonasymptotic convergence analysis for the unadjusted Langevin algorithm.  _Ann. Appl. Probab._ , 27(3):1551–1587, 2017\.
  * Ethier and Kurtz (1986) S. N. Ethier and T. G. Kurtz.  _Markov processes: Characterization and Convergence_.  Wiley Series in Probability and Mathematical Statistics. John Wiley & Sons, Inc., 1986.
  * Eyring (1935) H. Eyring.  The activated complex in chemical reactions.  _J. Chem. Phys._ , 3(2):107–115, 1935.
  * Fang et al. (1997) H. Fang, M. Qian, and G. Gong.  An improved annealing method and its large-time behavior.  _Stochastic Process. Appl._ , 71(1):55–74, 1997\.
  * Fournier and Tardif (2020) N. Fournier and C. Tardif.  On the simulated annealing in $\mathbb{R}^{d}$.  2020\.  arXiv:2003.06360.
  * Gao et al. (2020) X. Gao, Z. Q. Xu, and X. Y. Zhou.  State-dependent temperature control for Langevin diffusions.  2020\.  arXiv:2005.04507.
  * Ge et al. (2015) R. Ge, F. Huang, C. Jin, and Y. Yuan.  Escaping from saddle points – online stochastic gradient for tensor decomposition.  In _COLT_ , pages 797–842, 2015.
  * Geman and Hwang (1986) S. Geman and C.-R. Hwang.  Diffusions for global optimization.  _SIAM J. Control Optim._ , 24(5):1031–1043, 1986\.
  * Grenander and Miller (1994) U. Grenander and M. I. Miller.  Representations of knowledge in complex systems.  _J. Roy. Statist. Soc. Ser. B_ , 56(4):549–603, 1994.
  * Guo et al. (2020) X. Guo, J. Han, and W. Tang.  Perturbed gradient descent with occupation time.  2020\.  arXiv:2005.04507.
  * Holley et al. (1989) R. A. Holley, S. Kusuoka, and D. W. Stroock.  Asymptotics of the spectral gap with applications to the theory of simulated annealing.  _J. Funct. Anal._ , 83(2):333–347, 1989.
  * Jain and Kar (2017) P. Jain and P. Kar.  Non-convex optimization for machine learning.  _Found. Trends Mach. Learn._ , 10:142–336, 2017.
  * Jin et al. (2017) C. Jin, R. Ge, P. Netrapalli, S. Kakade, and M. I. Jordan.  How to escape saddle points efficiently.  In _ICML_ , pages 1724–1732, 2017.
  * Kirkpatrick et al. (1983) S. Kirkpatrick, J. Gelatt, and M. Vecchi.  Optimization by simulated annealing.  _Science_ , 220(4598):671–680, 1983.
  * Koulamas et al. (1994) C. Koulamas, S. Antony, and R. Jaen.  A survey of simulated annealing applications to operations research problems.  _Omega_ , 22(1):41–56, 1994.
  * Kramers (1940) H. A. Kramers.  Brownian motion in a field of force and the diffusion model of chemical reactions.  _Physica_ , 7(4):284–304, 1940.
  * Ma et al. (2019a) Y. Ma, N. Chatterji, X. Cheng, N. Flammarion, P. Bartlett, and M. I. Jordan.  Is there an analog of Nesterov acceleration for MCMC?  _arXiv:1902.00996_ , 2019a.
  * Ma et al. (2019b) Y. Ma, Y. Chen, C. Jin, N. Flammarion, and M. I. Jordan.  Sampling can be faster than optimization.  _Proc. Natl. Acad. Sci. USA_ , 116(42):20881–20885, 2019b.
  * Márquez (1997) D. Márquez.  Convergence rates for annealing diffusion processes.  _Ann. Appl. Probab._ , 7(4):1118–1139, 1997.
  * Menz and Schlichting (2014) G. Menz and A. Schlichting.  Poincaré and logarithmic Sobolev inequalities by decomposition of the energy landscape.  _Ann. Probab._ , 42(5):1809–1884, 2014.
  * Menz et al. (2018) G. Menz, A. Schlichting, W. Tang, and T. Wu.  Ergodicity of the infinite swapping algorithm at low temperature.  2018\.  arXiv:1811.10174.
  * Meyn and Tweedie (1993) S. P. Meyn and R. L. Tweedie.  Stability of Markovian processes. III. Foster-Lyapunov criteria for continuous-time processes.  _Adv. in Appl. Probab._ , 25(3):518–548, 1993\.
  * Miclo (1992) L. Miclo.  Recuit simulé sur $\mathbb{R}^{n}$. Étude de l’évolution de l’énergie libre.  _Annales de l’Institut Henri Poincaré_ , 28(2):235–266, 1992.
  * Miclo (1995) L. Miclo.  Une étude des algorithmes de recuit simulé sous-admissibles.  _Ann. Fac. Sci. Toulouse Math. (6)_ , 4(4):819–877, 1995.
  * Monmarché (2018) P. Monmarché.  Hypocoercivity in metastable settings and kinetic simulated annealing.  _Probability Theory and Related Fields_ , pages 1–34, 2018.
  * Otto and Villani (2000) F. Otto and C. Villani.  Generalization of an inequality by Talagrand and links with the logarithmic Sobolev inequality.  _J. Funct. Anal._ , 173(2):361–400, 2000.
  * Parisi (1981) G. Parisi.  Correlation functions and computer simulations.  _Nuclear Phys. B_ , 180(3, FS 2):378–384, 1981\.
  * Pavlyukevich (2007) I. Pavlyukevich.  Lévy flights, non-local search and simulated annealing.  _J. Comput. Phys._ , 226(2):1830–1844, 2007.
  * Raginsky et al. (2017) M. Raginsky, A. Rakhlin, and M. Telgarsky.  Non-convex learning via stochastic gradient Langevin dynamics: a nonasymptotic analysis.  In _COLT_ , pages 1674–1703, 2017.
  * Roberts and Tweedie (1996) G. O. Roberts and R. L. Tweedie.  Exponential convergence of Langevin distributions and their discrete approximations.  _Bernoulli_ , 2(4):341–363, 1996.
  * Rogers and Williams (1987) L. C. G. Rogers and D. Williams.  _Diffusions, Markov processes, and martingales. Vol. 2, Itô Calculus_.  Wiley Series in Probability and Mathematical Statistics. John Wiley & Sons, Inc., 1987.
  * Royer (2007) G. Royer.  _An initiation to logarithmic Sobolev inequalities_ , volume 14 of _SMF/AMS Texts and Monographs_.  American Mathematical Society, 2007.
  * Schlichting (2012) A. Schlichting.  _The Eyring-Kramers formula for Poincaré and logarithmic Sobolev inequalities_.  PhD thesis, Universität Leipzig, 2012.  Available at http://nbn-resolving.de/urn:nbn:de:bsz:15-qucosa-97965.
  * Shun and McCullagh (1995) Z. Shun and P. McCullagh.  Laplace approximation of high-dimensional integrals.  _J. Roy. Statist. Soc. Ser. B_ , 57(4):749–760, 1995.
  * Stroock and Varadhan (1979) D. W. Stroock and S. R. S. Varadhan.  _Multidimensional diffusion processes_ , volume 233 of _Grundlehren der Mathematischen Wissenschaften_.  Springer-Verlag, 1979.
  * Tang (2019) W. Tang.  Exponential ergodicity and convergence for generalized reflected Brownian motion.  _Queueing Syst._ , 92(1-2):83–101, 2019.
  * Tsybakov (2009) A. B. Tsybakov.  _Introduction to nonparametric estimation_.  Springer Series in Statistics. Springer, 2009.
  * van Laarhoven and Aarts (1987) P. J. M. van Laarhoven and E. H. L. Aarts.  _Simulated annealing: theory and applications_ , volume 37 of _Mathematics and its Applications_.  D. Reidel Publishing Co., 1987.
  * Vempala and Wibisono (2019) S. Vempala and A. Wibisono.  Rapid convergence of the Unadjusted Langevin Algorithm: isoperimetry suffices.  In _NeurIPS_ , volume 32, pages 8094–8106, 2019.
  * Wang et al. (2020) H. Wang, T. Zariphopoulou, and X. Y. Zhou.  Reinforcement learning in continuous time and space: A stochastic control approach.  _Journal of Machine Learning Research_ , 21:1–34, 2020\.
  * Zitt (2008) P. A. Zitt.  Annealing diffusions in a potential function with a slow growth.  _Stochastic Process. Appl._ , 118(1):76–119, 2008\.